Interferometric synthetic aperture microscopy

ABSTRACT

Methods and apparatus for three-dimensional imaging of a sample. A source is provided of a beam of substantially collimated light characterized by a temporally dependent spectrum. The beam is focused in a plane characterized by a fixed displacement along the propagation axis of the beam, and scattered light from the sample is superposed with a reference beam derived from the substantially collimated source onto a focal plane detector array to provide an interference signal. A forward scattering model is derived relating measured data to structure of an object to allow solution of an inverse scattering problem based upon the interference signal so that a three-dimensional structure of the sample may be inferred in near real time.

The present application claims the priority of U.S. Provisional PatentApplication Ser. No. 60/819,593, filed Jul. 10, 2006, which isincorporated herein by reference.

The present invention has been developed, in part, with Governmentsupport under Contract Numbers 02390265 and BES 05-19920, awarded by theNational Science Foundation, and Contract Number 1 R21 EB005321, awardedby the National Institutes of Health. The Government has certain rightsin the invention.

TECHNICAL FIELD

The present invention pertains to a computed imaging method forthree-dimensional interferometric determination of the susceptibility ofan object on the basis of coherently detected backscattered radiation,and, more particularly, to a method for providing improved resolutionoutside of the limited volume (the confocal volume) that can be resolvedin focus for typical optical, acoustic, and other microscopes.

BACKGROUND OF THE INVENTION

Optical coherence tomography (OCT), as employed for imaging throughscattering media, with significant applications for medical imaging, istypically based upon the notion of coherently ranging scatterers in afocal volume using a high-bandwidth, low temporal coherence source, suchas a superluminescent diode. By measuring the interferometriccross-correlation between a reference beam and the backscattered returnbeam, the scatterers corresponding to a particular time delay, typicallycorresponding to a depth in the medium, can be ranged. Techniquesemploying OCT are surveyed in Bouma et al., Handbook of OpticalCoherence Tomography (Marcel Dekker, 2001), which is incorporated hereinby reference. Optical coherence microscopy (OCM), described, forexample, by Schmitt et al., Subsurface Imaging of Living Skin withOptical Coherence Microscopy, Dermatology, vol. 191, pp. 93-98 (1995)and Izatt, Optical Coherence Microscopy in Scattering Media, OpticsLetters, vol. 19, pp. 590-592 (1994), incorporated by reference herein,combines the depth-ranging capability of OCT with confocal microscopy toobtain micron-scale imaging beneath the surface of a scattering medium.

Traditional OCT operates by illuminating the object with a focused,broadband beam. The back-scattered light is collected, and by usinginterferometric detection, the time delay and therefore the distancealong the beam to scatterers inside the object is determined. Byscanning the beam through the object, the locations of scatterers inthree dimensions can be found. The resolution of OCT in the axialdirection, along the direction of propagation of the beam, is determinedprimarily by the bandwidth of the light source. However, the resolutionin the transverse direction is not constant along the beam. At the focusof the beam, the resolution is determined by the focused spot size, butaway from the focus, the resolution degrades because the beam isdiverging or converging. This loss of resolution is usually assumed tobe an inevitable consequence of defocus. Because the depth-of-field orthe confocal volume decreases in size as the numerical aperture of theillumination and detection optics increase, improved transverseresolution results, according to current practice, in a smaller range ofdepths that can be resolved for a typical OCT system unless the focus ismechanically scanned. “Confocal volume,” as used herein, refers to avolume that is essentially in focus, within a specified criterion, at asingle relative placement of the sample relative to the optical system.

SUMMARY OF THE INVENTION

In accordance with preferred embodiments of the present invention, amethod is provided for determining the three-dimensional susceptibilityof a sample. The method has steps of:

-   -   a. providing a source of a beam of radiation characterized by a        temporally dependent spectrum and a local propagation axis;    -   b. irradiating the surface in a plane characterized by a fixed        displacement along the propagation axis of the beam;    -   c. superposing scattered radiation from the sample with a        reference beam derived from the source of the beam to provide an        interference signal;    -   d. deriving a forward scattering model relating measured data to        structure of an object; and    -   e. solving an inverse scattering problem based upon the        interference signal to infer a three-dimensional susceptibility        of the sample.

In accordance with alternate embodiments of the invention, the forwardscattering model may relate a transform of measured data to a transformof structure of the sample. A step may additionally be provided oftransforming the interference signal into a transform of theinterference signal in a transform domain. The step of solving aninverse scattering problem may be performed in substantially real time.

In accordance with further embodiments of the invention, the step ofirradiating may include focusing radiation substantially at the surfaceof the sample, such as through a microscope objective. It may alsoinclude delivering the beam to the irradiated sample by at least one ofa catheter, a needle, a probe, an endoscope, or a laparoscope. Themethod may also entail sweeping the spectrum of the source of radiationas a function of time.

In other embodiments of the invention, a tunable laser may be providedas the source of radiation. The step of superposing scattered radiationfrom the sample with a reference beam may include configuring arms ofthe beam in an interferometer, such as a Michelson a Mach-Zehnder, and aFizeau interferometer.

In accordance with another aspect of the invention, an interferometricsynthetic aperture microscope (ISAM) is provided that has a source ofsubstantially coherent illumination for illuminating an object subjectat a fixed focal distance relative to a fiducial position defined withrespect to the microscope. The ISAM also has a reference arm forrelaying a portion of the coherent illumination onto a focal planesensor, a telescope for relaying a field scattered by an object onto thefocal plane sensor, and a processor for deconvolving the field measuredby the focal plane sensor to recover a three-dimensional image of theobject. The reference arm, in particular, may be a fixed referencemirror.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing features of the invention will be more readily understoodby reference to the following detailed description, taken with referenceto the accompanying drawings, in which:

FIG. 1 is an exemplary schematic depiction of an optical coherencetomography system in its time-domain and spectral-domain variants.

FIG. 2( a) is an example of OCT data derived from microparticlesembedded in a tissue phantom outside of the focus, while FIG. 2( b) isan ISAM reconstruction of the data.

FIG. 3( a) is a schematic depiction of a full-field OCT system inaccordance with a preferred embodiment of the present invention.

FIG. 3( b) depicts an example of a band volume for an instrument, inaccordance with the present invention, having a numerical aperture of0.5 and bandwidth from 0.8 k_(max)<k<k_(max).

FIG. 3( c) shows the projection of the simulated data collapsed (summed)along one transverse direction, while FIG. 3( d) is a projection of thecomputed reconstruction of the scatterers.

FIGS. 4( a)-4(c) show three pairs of en face images of the time-domaindata (left side) and the reconstructed volume (right).

FIG. 5( a) compares a volume isosurface plot of the raw data with thereconstructed computed image of FIG. 5( b).

FIG. 6 is a schematic depiction of an interferometric synthetic aperturemicroscope using spectral detection in accordance with an embodiment ofthe present invention.

FIG. 7 is a schematic view of an OCT catheter shown in cross-section.

FIG. 8 defines the coordinates used in the description of the OCTcatheter.

FIG. 9( a) is simulated OCT for randomly scattered point objects andFIG. 9( b) is a reconstruction of the point sources from the simulateddata.

FIG. 10 is a plot depicting full-width-half-maximum transversepoint-spread-function (PSF) resolution of simulated point sourcessituated at different distances from the catheter axis, as a function offocus radius and beam width. The abscissa is the distance from theorigin from which the simulated point source is placed. In each part,the number next to each curve is the distance away from the origin thebeam was focused. The beam waist for each part is (a) 1.5λ,(b) 3λ,(c)4.5λ, and (d) 6λ.

FIG. 11 is a plot depicting full-width-half-maximum axial resolution ofsimulated point sources imaged with two different beam widths focused 25wavelengths from the catheter axis, for various fractional bandwidths ofthe source. The dotted curve corresponds to a 1.5λ beam waist, while thesolid curve corresponds to a 5λ beam waist.

FIG. 12 shows interferometric data from a tissue phantom consisting oftitanium dioxide scatterers suspended in silicone. Planar slices the 3-Dunprocessed data (left) and ISAM reconstruction (right) are shown fortwo en face planes above the focus and one below the focus. Magnifiedunprocessed sections for three depths are shown in (a) at z=1100 μm, (b)at z=475 μm, and (c) at z=−240 μm, where z=0 μm is the focal plane.Magnified ISAM reconstructions for these corresponding planes are shownin (d), (e), and (f), respectively. The scale bar represents 18 μm.

FIGS. 13( a)-(c) show en face images from human breast tissue. (a)Unprocessed interferometric data and (b) ISAM reconstructions are shownfor depths located at z=591 μm (Slice A) and z=643 μm (Slice B), wherez=0 μm is the focal plane. (c) The corresponding histological sectionsare shown for comparison.

FIG. 14 is a flowchart depicting steps in imaging a sample using inaccordance with embodiments of the present invention.

FIG. 15( a) depicts the relation between spatial frequencies of thesignal space and spatial frequencies in the object space; FIG. 15( b)shows a sampling lattice for selected (β, |Q|) values on a uniform (k,|Q|) grid; and FIG. 15( c) shows a sampling lattice for selected (k,|Q|) values on a uniform (β, |Q|) grid.

FIG. 16 is a data flow chart for ISAM processing, in accordance withembodiments of the present invention.

FIG. 17 is a computational flow chart for memory allocation forsuccessive steps of ISAM processing, in accordance with embodiments ofthe present invention.

DETAILED DESCRIPTION OF SPECIFIC EMBODIMENTS

The following detailed description may be more fully understood byreference to Ralston et al., Inverse Scattering for Optical CoherenceTomography, J. Opt. Soc. Am. A, vol. 23, pp. 1027-37 (May, 2006), whichis appended hereto, and incorporated herein by reference.

FIG. 1 shows an example of an OCT system, designated generally bynumeral 5, including options to detect the cross-correlation signal ineither the spectral domain or the time domain. This particular OCTsystem uses a coherent light source 12, such as a neodymium:vanadate(Nd:YVO₄) pumped titanium:sapphire laser. In the exemplary embodimentdepicted, a center wavelength of 800 nm and an average bandwidth of 100nm yields an axial resolution of ˜2 μm in tissue. The beam 14 islaunched, via coupling optics 15, into a fiber-optic cable 16 andconnected to a pair of fiber-optic beam splitters 18, 20 (Gould FiberOptics, Inc., Millersville, Md.). The interferometer employs asingle-mode 50/50 fiber optic splitter 20 that delivers and couples thesignals to and from the sample arm 22 and the galvanometer-basedreference arm 24. The 90/10 beam splitter 18 serves the purpose ofproducing a reference spectrum for source noise reduction in thebalanced photodetector. This setup can make use of both the singlephotodetector and the high-speed dual balanced photodetectors 26 and 28.In the sample arm, interchangeable achromatic lenses 30, from 9 mm to 40mm, can be used to focus from 10 mW to 30 mW of light down tocorresponding spot size (determined by transverse resolution).Polarization paddles 32 in the sample arm are able to change thepolarization incident upon the sample to achieve maximum interference.

A time-domain system is characterized by capturing axial scans while areference arm 24 varies the path length in relation to the sample arm22. There are several designs of varying the path length in thereference arm. In certain embodiments, reference arms may include astandard spatial domain delay, or a rapid scanning optical delay (RSOD)which employs a spatial Fourier shift to create a delay. Typically, asingle or dual-balanced photodetector 36 is used to capture interferencedata.

A spectral-domain system is characterized by capturing the spectrum of alight source and inferring the spatial axial susceptibility.Susceptibility describes, in theoretical terms, the dielectricproperties of a material that give rise to scattering and thus to therepresentation ordinarily considered to constitute an image. There areseveral ways to detect the spectrum including employing a frequencyswept laser and detection with a single photodetector, or a broadbandlaser and detection using an imaging grating spectrometer 40 so that allfrequencies are detected simultaneously.

In order to elucidate the present invention, a theory of inversescattering is presented that has been developed for optical coherencetomography (OCT) and that is used to resolve three-dimensional objectstructure, taking into account the finite beam width and focusing. Whilethe invention is described herein in optical terms, it is to beunderstood that the invention is not so limited, and that the teachingsprovided herein are equally applicable to any radiation, whetheracoustic or of other particles, massive or otherwise, that may becharacterized in terms of wave propagation.

By using the invention, transverse and axial resolution produced byconventional OCT imaging inside the confocal volume can be achievedoutside of the confocal volume. Explicitly, experiments show thatscatterers can be resolved outside of the confocal volume withresolution comparable to that achievable inside the confocal volume.Numerical simulations and experimental results demonstrate theeffectiveness of this technique. When the algorithm is applied toexperimentally-acquired OCT data, the transverse and axial resolutionsoutside of the confocal parameter are improved, extending the apparentconfocal parameter range. These results further improve thehigh-resolution cross-sectional imaging capabilities of OCT.

To illustrate the problem in OCT that is solved with ISAM, a sample isdesigned and imaged, which clearly shows the effect of the probing beam.A tissue phantom, a collection of titanium dioxide scatterers having amean diameter of 1 μm suspended in silicone, was imaged with aspectral-domain OCT (SD-OCT) system. FIG. 2( a) displays the OCT datawith no beam consideration. The 2D ISAM reconstructed image of an areaof 500 pm (transverse) by 1000 μm (axial) is depicted in FIG. 2( b),where the bandwidth is 100 nm, the focal length of the lens is 15 mm,the spot size is 11 μm, and the confocal parameter is 996 μm. The imageresolution of point scatterers outside of the confocal region for theoriginal experimental image data is not constant, but for thereconstruction, the resolution is relatively constant along the entireimage with only amplitude variations. The interference between the lightscattered from a group of adjacent particles (boxed) is evident in theoriginal image (top magnified). Methods in accordance with the presentinvention properly rephase the signal from scatterers to produce awell-resolved image (bottom magnified).

The sample imaged in FIG. 2 consists of a uniform number of scatterers,yet when looking at the reconstruction the number of scatters outside ofthe confocal region seems to increase. The 2D equivalent of thealgorithm resolves only a perpendicular projection of the 3D data set.The imaging beam has a wider contour outside of the confocal region, andthus the beam is incident upon a larger volume illuminating morescatterers.

ISAM fundamentally relies on the solution of the inverse scatteringproblem, S=Kη, where K is a linear operator which transforms thecollected signal S to the object's spatial susceptibility η. (Thedimensionalities of S and η are discussed below.) In some geometries, Kcan be diagonalized in the Fourier domain producing an efficientalgorithm for attaining spatially invariant resolution.

There are many embodiments for ISAM, each of which relies on specificimaging geometries and the various embodiments may advantageously beemployed in variations on the instrumentation. Different geometriesrequire somewhat different formulations of the forward and inversescattering problems. Included here are the solutions for threegeometries that are most likely to be useful in practice. The firstgeometry, for a Gaussian beam with a focus transversely scanned over aplane, is the most often used geometry of OCT when external surfaces areto be imaged. The second geometry, for an azimuthally scanned beam, istypically employed for catheters for imaging internally inside the humanbody, and may be utilized in applications wherein the beam is deliveredto the irradiated sample by a needle, a probe, an endoscope, or alaparoscope. Finally, the full-field geometry is useful when externalsurfaces are to be imaged and speed of data acquisition is paramount.Other instruments will likely be amenable to be adapted to thesegeometries, and the ability to perform ISAM is not limited to thesegeometries. The types of incident fields that we have specificallyaddressed are focused Gaussian beams and plane waves because they arethe most common and readily produced and detected beams used for OCT,however the scope of the present invention is not limited to these typesof incident field.

In accordance with embodiments of the present invention, thecapabilities of both OCT and OCM are greatly extended by computedimaging and synthetic aperture techniques. Among recently demonstratedadvantages is the ability to resolve features in the sample that areoutside of the confocal region. A more quantitatively accurate andfaithful representation of the sample structure than hitherto availableis provided by solution of the inverse scattering problem as appliedboth to full-field planar OCT/OCM as well as to OCT from anazimuthally-scanned catheter. In either case, and in accordance withpreferred embodiments of the invention, the focus may advantageouslyremain fixed at the surface of the sample, while computed imagingtechniques are used to infer the structure inside and outside of thedepth-of-field. By applying the invention, the focus need not be scannedthrough the sample.

As further described below, a forward scattering model is derived whichrelates the measured data to the object structure. From this model, asolution of the inverse scattering problem is obtained that infers theobject structure from the data. The achievable resolution and systembandpass is also derived from this forward model, and application of themethod is demonstrated by means of a simulation.

Full-Field Non-Scanned Beam Implementation

By means of the novel methods described herein, computed imagingtechniques are employed to reconstruct features that are both inside andoutside the focus. Instead of scanning the focus through the sample, thefocus is fixed at the surface of the sample, and no relative translationis needed between the objective and the sample. A frequency-swept sourcecan be utilized to provide a new degree of freedom, replacinginformation lost by fixing the focus, without sacrificing detail outsideof the focus. Because the objective and sample can remain fixed relativeto each other, no translation hardware is needed which makes placing theobjective on a fiber optic or a handheld probe easier. This method maylead to faster and more accurate full-field OCT imaging because dataacquisition can be very rapid, requiring only that the two-dimensionalinterferogram is sampled while the frequency of the source is scanned.By virtue of computational image formation, the need to physically forman image of each plane in the volume, as typically done in full-fieldOCT, is obviated. As data acquisition speed and computational speedcontinue to increase, video-rate scanning of three-dimensional volumesmay become possible.

In order to provide an understanding of computational image formation inthe context of full-field OCT, a physical model for the scatteringprocess is developed, and from this, a relationship between the data andthe object structure is derived. Based on this relationship, the inversescattering problem is solved to infer the sample structure from thedata.

Full-field OCT allows an entire plane of scatterers to be rangedsimultaneously, which makes it a very rapid way to acquire the structureof a volume. A full-field OCT system that is typical of the currentstate-of-the-art consists of a Michelson interferometer, again, with abroadband illumination source. Reference and sample beams are derivedfrom the broadband source using a beam splitter. An extended area of thesample is illuminated by a broadband collimated beam through amicroscope objective. The objective is focused at the depth of featuresof interest. A signal is scattered by the sample back through theobjective. A reference beam is delayed to return to the beam splitter atthe same time the signal scattered from the focus arrives. The referenceand signal are superimposed and focused on a focal plane array (such asa charge-coupled device (CCD) sensor) which detects the interferencesignal. The amplitude of the interference signal corresponds to thereflections of scatterers at the focus plane. By translating the samplethrough the focus plane, the scatterers at many different depths may beranged.

While this method can be used to obtain high resolution images forentire volumes of a sample quickly, it has a number of disadvantages.First, the sample and microscope objective must be translated relativeto each other, which is relatively slow and requires fine positioning.Second, this method uses time-domain detection that produces a lowersignal-to-noise ratio than Fourier-domain, or frequency-swept, OCT.

A full-field OCT system, in accordance with embodiments of the presentinvention, is now described with reference to FIG. 3( a) and isdesignated generally by numeral 100. While system 100, as shown, isbased on a Michelson interferometer, other interferometricconfigurations such as that of a self-referencing Fizeau design, may beused and are encompassed within the scope of the present invention andof any claims appended hereto. In system 100, the illumination source isa tunable narrow band laser 112. It is to be understood that partiallycoherent sources may also be employed within the scope of the presentinvention, where their application is consistent with the methodsdescribed, and that references herein to a laser may also encompasssources that lack temporal or spatial coherence, or both, unless theparticular context dictates otherwise.

Laser 112 is tuned to wavelengths λ that correspond to spatialfrequencies k. Laser 112 nominally emits a plane wave (or is spatiallyfiltered to produce one). The coherence length of this laser should beat least as twice as long as the total depth of the sample 8 understudy, to ensure that fields scattered throughout the entire samplesimultaneously interfere with the reference field.

Laser illumination 113 is split by a beam splitter 114 into twocomponents. One component 115 travels to a reference mirror (or delaymirror) 116, and is reflected back through the beamsplitter 114 to theoutput port where the focal plane array 108 is located. It is to beunderstood that, most generally, numeral 108 designates a detector, andthat references to detector 108 as a focal plane array are by way ofnon-limiting example. The other beam component 117 is demagnified by afactor 1/M using a telescope 118 of magnification M. The purpose oftelescope 118 is to concentrate the illumination onto the sample 8, andthen relay a magnified scattered field 120 to the focal plane array 108.This telescope consists of two converging lenses: a relay lens 122 and amicroscope objective 124. The illumination on the sample is a normallyincident plane wave 126. The sample scatters some radiation 128backwards through the telescope 118. The telescope is aligned toafocally and telecentrically image the front surface of the sample tothe focal plane array. Sitter et al., Three-dimensional Imaging: aSpace-invariant Model for Space Variant Systems, Appl. Opt., vol. 29,pp. 3789-94 (1990) discusses three-dimensional imaging problems, and isincorporated herein by reference.

It is to be noted, significantly, that in a manner distinct from that ofstandard full-field OCT microscopy, the focus of the objective 124 mayremain fixed, in accordance with the present invention, at the surfaceof sample 8. For purposes of the following heuristic analysis, it isassumed that telescope 118 is aberration free and vignetting inside thetelescope is negligible. If the telescope is assumed to correctspherical aberration, then there is a finite volume within the samplespace for which these assumptions hold. A pupil 130 is placed at thefocus of the illumination beam inside the telescope to spatially filterthe backscattered signal so as to enforce a well-defined spatialbandlimit. The response of the telescope is modeled by a space-invariantconvolution with a bandwidth determined by the pupil size. At the focalplane array 108, the reference and sample signals superimpose andinterfere, and the intensity of the interference is detected. Theintereference signal from detector 108 is coupled to an InverseScattering Solver 132, the operation of which is now described.

To derive the relationship between the object structure and the datadetected on the sensor, a mathematical model of scattering of theillumination field by the object and interferometric detection at thesensor is now described with reference to FIG. 14. For convenience ofdescription, a scalar field is substituted for the electromagneticfield, neglecting polarization effects. The incident field on thesample, provided in step 140, is given by the expression:E _(i)(r;k)=A(k) exp(ikz)  (1)where r is a location within the volume of sample 8, k is the spatialfrequency of the illumination, A(k) is the amplitude of the illuminationat frequency k, and {circumflex over (z)} is the direction of increasingdepth into the sample. For present purposes, it is assumed that thescattering is well-modeled by the weak or first Born scatteringapproximation, where the scattering by the object is modeled as asource. The susceptibility of the object is given by η(r) such thatη(r)=0 for z<0.

The secondary scattered field E_(s)(r′; k) from the object at the planez=0 is given by the expression

$\begin{matrix}{{E_{s}\left( {r^{\prime};k} \right)} = {\int\limits_{V}{{\mathbb{d}^{3}{{rE}_{i}\left( {r;k} \right)}}{\eta(r)}{\frac{\exp\left( {{\mathbb{i}}\; k{{r^{\prime} - r}}} \right)}{{r^{\prime} - r}}.}}}} & (2)\end{matrix}$

To simplify this relationship, it is useful to define thetwo-dimensional Fourier transform

${{\overset{\sim}{E}}_{s}\left( {q;k} \right)} = {\int\limits_{z^{\prime} = 0}{{\mathbb{d}^{2}r^{\prime}}{E_{s}\left( {r^{\prime};k} \right)}{\exp\left( {{\mathbb{i}}\;{q \cdot r^{\prime}}} \right)}}}$with q being a transverse spatial frequency such that q·{circumflex over(z)}=0. In terms of q, Eq. (2) is found to be

$\begin{matrix}{{{\overset{\sim}{E}}_{s}\left( {q;k} \right)} = {{A(k)}{\int_{V}{{\mathbb{d}^{3}r}\;{\eta(r)}\exp\left\{ {{{\mathbb{i}}\left\lbrack {q \cdot r} \right\rbrack} + {{\mathbb{i}}\;{z\left\lbrack {k + {k_{z}(q)}} \right\rbrack}}} \right\}{k_{z}(q)}^{- 1}}}}} & (3)\end{matrix}$where k_(z)(q)=√{square root over (k²−q²)}, and substituting Eq. (1)into Eq. (2). (In accordance with convention, x² designates the squaremodulus of a vector x.) The three-dimensional Fourier transform isdefined such that

${\overset{\sim}{\eta}(Q)} = {\int\limits_{V}{{\mathbb{d}^{3}r}\;{\eta(r)}{{\exp\left( {{\mathbb{i}}\;{Q \cdot r}} \right)}.}}}$It is then found that the right-hand integral can be expressed in termsof {tilde over (η)}(Q):{tilde over (E)} _(s)(q;k)=A(k)k _(z)(q)⁻¹ {tilde over(η)}{q+{circumflex over (z)}[k+k _(z)(q)]}  (4)

The field E_(f)(r;k) is produced by the propagation of E_(s)(r′;k)through telescope 118 to focal plane array 108 (shown in FIG. 3( a)).Because the telescope is assumed to be an aberration-free telescopewhich afocally and telecentrically images the plane at the sample z=0 tothe focal plane array in the plane z=z_(f), its function can be modeledby a simple convolution with a point spread function accounting for thefinite bandwidth of the telescope, and a magnification factor given byM. The field at the focal plane array is given by E_(f)(r;k), and thepoint spread function of the telescope is given by P(r;k). Therelationship between E_(f)(r;k) and E_(s)(r′;k) isE _(f)(r;k)=M ⁻¹ ∫d ² r′E _(s)(r′;k)P(r/M−r′;k)  (5)where the factor M⁻¹ accounts for energy conservation between the twoplanes.

${{\overset{\sim}{E}}_{f}\left( {q;k} \right)} = {\int\limits_{z = z_{f}}{{\mathbb{d}^{2}{{rE}_{f}\left( {r;k} \right)}}{\exp\left( {{\mathbb{i}}\;{q \cdot r}} \right)}}}$and the coherent transfer function of the telescope {tilde over(P)}(q;k)=∫d²r P(r;k)exp(iq·r), the convolution of Eq. (5) may beexpressed as{tilde over (E)} _(f)(q;k)=M{tilde over (E)} _(s)(Mq;k){tilde over(P)}(Mq;k)=MA(k){tilde over (P)}(Mq;k)k _(z)(Mq)⁻¹ {tilde over(η)}{Mq+{circumflex over (z)}[k+k _(z)(Mq)]}  (6)Eq. (6) specifies a relationship between Fourier components of the fieldon the focal plane array and those of the object susceptibility.

In accordance with preferred embodiments of the invention, referencemirror 116 is placed to effect a delay on the reference beam 145 of τrelative to the total delay required for the beam to travel from thebeamsplitter 114 to the plane z=0 in the sample arm and back. Thereference field E_(r) (r; k, τ), relayed to the focal plane array isthen given byE _(r)(r;k,τ)=A(k)exp[iω(k)τ],  (7)where ω(k) is a dispersion relation relating the temporal frequency withthe spatial frequency in the sample medium.

For example, if the object medium is free space, then ω(k)=kc, where cis the speed of light in vacuum. The interference intensityI(r;k,τ)=|E_(r)(r;k,τ)+E_(f)(r;k)|² on the focal plane array is thengiven byI(r;k,τ)=|A(k)² +|E _(f)(r;k)²+2A(k)Re{E _(f)(r;k)exp[−iω(k)τ]}.  (8)The part of the signal 142 that is due to interference between thesignal and reference beams occurring in step 141 is defined as the datafunction D(r;k)=A(k)E_(f)(r;k). The complex D(r;k) can be estimated frommeasurements of I(r;k,τ). For example, three measurements of I(r;k,τ)such that ωτ=0, π/2, and π may be summed (in step 147 of FIG. 14) toyield, as an approximation:

$\begin{matrix}{{D\left( {r;k} \right)} = {{\frac{1 - {\mathbb{i}}}{4}{I\left( {{r;k},0} \right)}} - {\frac{1 + {\mathbb{i}}}{4}{I\left( {{r;k},{\pi/\omega}} \right)}} + {\frac{\mathbb{i}}{2}{{I\left( {{r;k},{{\pi/2}\omega}} \right)}.}}}} & (9)\end{matrix}$

The foregoing method of phase-shifting for interferometry is described,for example, in Hariharan, Optical Interferometry (Academic Press,2003), which is incorporated herein by reference. Inserting the resultsof Eq. (6), the Fourier transform of the data function, which is {tildeover (D)}(q;k)=∫d²r D(r;k)exp(iq·r), can be expressed as{tilde over (D)}(q;k)={tilde over (K)}(q;k){tilde over(η)}{Mq+{circumflex over (z)}[k+k _(z)(Mq)]},  (10)where, for convenience, the bandpass function is defined{tilde over (K)}(q;k)=MA(k)² {tilde over (P)}(Mq;k)k _(z)(Mq)⁻¹  (11)Thus, the data are expressed in terms of the 3-D Fourier transform ofthe sample structure, and, so, the resolution of the reconstruction ofthe sample structure is space invariant. However, vignetting andaberrations in the telescope limit the volume over which this resolutioncan be obtained.

To obtain the measurements needed to reconstruct η(r), one must varyboth k and τ. In practice, however, it is often slow and inconvenient toadjust both. If one is willing to tolerate some image artifacts, justone of these parameters need be scanned. For simplicity, it is assumedthat the pupil function P(r′;k) is real and symmetric, which is oftenthe case (for example, when pupil 130 is circular), so that {tilde over(P)}(q;k) is likewise real and symmetric.

If mirror 116 is fixed such that τ=0, then the imaginary component ofD(r;k) is not obtainable. If the imaginary part of D(r;k) is assumed tobe zero, then due to the Hermitian symmetry of the Fourier transform ofreal functions {tilde over (D)}(−q,k)={tilde over (D)}(q,k)*. Thefunction {tilde over (η)}(Q) then also has Hermitian symmetry reflectedover the axis. The effect is that a conjugate image of thesusceptibility is present, reflected across the plane z=0. Because thedelay τ=0 corresponds to the plane z=0, as long as the entire sample iscontained such that z>0, the conjugate image and the real image do notoverlap. In addition, there is an artifact corresponding to the term|E_(f)(r;k)|² in Eq. (8). If the magnitude of the sample signal is smallrelative to the reference signal, the magnitude of this artifact is alsosmall compared to the real image and can be neglected. This term mayalso be eliminated by modulating the phase of the reference field andlocking in only on the modulation, i.e., by phase-sensitive detection ofthe intereference signal.

If the delay τ is scanned, and the laser emits all frequencies ksimultaneously (such as occurs in a mode locked laser or a spontaneousemission source), the signal I_(T)(r;τ) is the sum of the interferencepatterns over all emitted frequencies:

$\begin{matrix}{{I_{T}\left( {r;\tau} \right)} = {{\frac{1}{2\;\pi}\left\lbrack {\int_{- \infty}^{\infty}\ {{\mathbb{d}{k\left( \frac{\mathbb{d}\omega}{\mathbb{d}k} \right)}}\left( {{{A(k)}}^{2} + {{E_{f}\left( {r;k} \right)}}^{2}} \right)}} \right\rbrack} + {\frac{1}{\pi}{Re}{\left\{ {\int_{- \infty}^{\infty}\ {{\mathbb{d}{k\left( \frac{\mathbb{d}\omega}{\mathbb{d}k} \right)}}{D\left( {r;k} \right)}{\exp\left\lbrack {{- {\mathbb{i}}}\;{\omega(k)}\tau} \right\rbrack}}} \right\}.}}}} & (12)\end{matrix}$

The term in square brackets in Eq. (12) is a background intensity thatis independent of τ and therefore is easily subtracted to remove itscontribution from the measured intensity. Neglecting the backgroundintensity and the slowly-varying Jacobian

$\left( \frac{\mathbb{d}\omega}{\mathbb{d}k} \right),$Eq. (12) relates the real part of the inverse Fourier transform ofD(r;k) with respect to k to the total intensity I_(T)(r;τ). To be ableto remove the Re{ } operation so that a unique solution for D(r;k) canbe found, one equates D(r;−k)=D(r;k)*. Eq. (10) then likewise enforcesHermitian symmetry on η(−Q)=η(Q)*. Therefore in this case thereconstructed susceptibility is assumed to be real-valued.

In this derivation, the focal plane of the objective and the frontsurface of the sample are assumed to coincide (at z=0). This assumptionhas simplified the preceding analysis and presentation, but it is notrequired within the scope of the present invention. If the focus isplaced below the sample surface by a distance z₀, but the delay producedby the reference still coincides with the delay of the sample surface,the data can be modified to account for the displacement. In particular,the modified data {tilde over (D)}′(q; k) is related to the sampled data{tilde over (D)}(q; k) by:{tilde over (D)}′(q;k)={tilde over (D)}(q;k)exp{iz ₀ [k−k_(z)(Mq)]}  (13)This formula can be found by noting that the field relayed by thetelescope is now situated at the plane z=z₀, adding an extra termexp{−iz₀[k+k_(z)(q)]} to the right side of Eq. (3). At the same time,the delay reference mirror must be moved a distance further from thebeamsplitter so that the new effective delay corresponds to the frontsurface of the sample, adding a factor of exp(−2ikz₀) to the right sideof Eq. (7) to place the reference delay coincident with the frontsurface of the sample. Effectively the measured field is computationallypropagated at each frequency to the surface of the sample.

Using the mathematical model 146 developed in the foregoing discussion,a solution to the inverse scattering problem may be derived in step 143.In general, the solution is ill-posed and so regularization techniquesare used to produce a stable solution. Because the forward problem islinear, we derive a linearized inverse based on least-squares error. Todo so, we first specify the complete forward operator K such that D=Kη,which relates the data to the object structure{tilde over (D)}(r;k)=Kη=∫d ³ r′K(r′,r;k)η(r′)  (14)where the kernel K(r′,r;k) of the operator K is given by

$\begin{matrix}{{K\left( {r^{\prime},{r;k}} \right)} = {M^{- 1}{A(k)}^{2}{\int{{\mathbb{d}^{2}r^{''}}\frac{\exp\left( {{\mathbb{i}}\; k{{r^{''} - r^{\prime}}}} \right)}{{r^{''} - r^{\prime}}}{{P\left( {{{r/M} - r^{''}};k} \right)}.}}}}} & (15)\end{matrix}$

Given this relation between the data and the object, we can define aleast-squares solution η⁺ (r) for object susceptibility as

$\begin{matrix}\begin{matrix}{{\eta^{+}(r)} = {\underset{\eta}{\arg\;\min}{{D - {K\;\eta}}}^{2}}} \\{= {\underset{\eta}{\arg\;\min}{\int{{\mathbb{d}^{2}r^{\prime}}{\int{{\mathbb{d}k}{{{{\overset{\sim}{D}\left( {r^{\prime};k} \right)} - {K\;{\eta(r)}}}}^{2}.}}}}}}}\end{matrix} & (16)\end{matrix}$

Expressed in operator notation, the solution to this least squaresproblem is given by the pseudo inverse η⁺=(K^(†)K)⁻¹K^(†)D where K^(†)is the Hermitian conjugate of K and K^(†)K is assumed to be invertible.It is much simpler to formulate the least-squares problem in the Fourierdomain, using the relation of Eqs. (10) and (11). If we instead definethe operator K such that D=K{tilde over (η)}. This operator can be usedto construct a least squares solution {tilde over (η)}⁺ such that:

$\begin{matrix}\begin{matrix}{{{\overset{\sim}{\eta}}^{+}(Q)} = {\underset{\overset{\sim}{\eta}}{\arg\;\min}\left( {{{D - {K\;\overset{\sim}{\eta}}}}^{2} + {\gamma{\overset{\sim}{\eta}}^{2}}} \right)}} \\{= {\underset{\overset{\sim}{\eta}}{\arg\;\min}\left( {\int{{\mathbb{d}^{2}q}{\int{{\mathbb{d}k}{{{\overset{\sim}{D}\left( {q;k} \right)} - {{\overset{\sim}{K}\left( {q;k} \right)}\overset{\sim}{\eta}}}}}}}} \right.}} \\\left. {{\left\{ {{Mq} + {\hat{z}\left\lbrack {k + {k_{z}({Mq})}} \right\rbrack}} \right\} }^{2} + {\gamma{\overset{\sim}{\eta}}\left\{ {{Mq} + {\hat{z}\left\lbrack {k + {k_{z}({Mq})}} \right\rbrack}} \right\}^{2}}} \right)\end{matrix} & (17)\end{matrix}$with {tilde over (K)}(q;k) taken from Eq. (11). A Tikhonovregularization term with regularization constant γ has been added tostabilize the solution and ensure that a unique solution exists. Thesolution {tilde over (η)}^(†) is given in step 144 by

$\begin{matrix}\begin{matrix}{{{\overset{\sim}{\eta}}^{+}\left\{ {{Mq} + {\hat{z}\left\lbrack {k + {k_{z}({Mq})}} \right\rbrack}} \right\}} = {\left( {{K^{\dagger}K} + \gamma} \right)^{- 1}K^{\dagger}D}} \\{= {\frac{{\overset{\sim}{D}\left( {q;k} \right)}{{\overset{\sim}{K}}^{*}\left( {q;k} \right)}}{{{\overset{\sim}{K}\left( {q;k} \right)}}^{2} + \gamma}.}}\end{matrix} & (18)\end{matrix}$Resolution and Bandpass of Full-Field ISAM

Eq. (10) expresses a relationship between the 2-D Fourier transform ofthe data and the 3-D Fourier transform of the object. As mentionedpreviously, this relationship implies that the resolution of thereconstructed object is space invariant. With suitable specifications ofthe instrument, it is possible to identify the region of the Fourierspace of the structure function that can be sampled. This region iscalled the “band volume” and is an analogue to the bandlimit ofone-dimensional signals, except that the band volume consists of theinterior of a shape in three-dimensional Fourier space rather than justa one-dimensional interval.

There are two specifications of the instrument that determine the shapeof the band volume. The first is the bandwidth of the illumination,which is specified by the interval of frequencies k_(min)<k<k_(max). Theother parameter is the numerical aperture (NA) of the imaging system0<NA<1. A particular numerical aperture implies a pupil bandpass{tilde over (P)}(q;k)=1 for|q|≦(NA)k{tilde over (P)}(q;k)=0 for|q|>(NA)k.  (19)

These inequalities constrain which points on the data function {tildeover (D)}(q;k) can be sampled. The relation of Eq. (10) is a one-to-onemapping between each of these points in the data function and points inthe 3-D Fourier space of the object {tilde over (η)}(Q). The band volumeis the intersection of the volumes defined by the two inequalitiesexpressed in terms of the object spatial frequency Q

$\begin{matrix}{k_{\min} < {Q^{2}/\left( {2{Q \cdot \hat{z}}} \right)} < {{k_{\max}\left( {2{Q \cdot \hat{z}}} \right)}{\sqrt{Q^{2} - \left( {Q \cdot \hat{z}} \right)^{2}}/Q^{2}}} < {{NA}.}} & (20)\end{matrix}$

FIG. 3( b) shows an example of a band volume for an instrument with 0.5NA and bandwidth from 0.8 k_(max)<k<k_(max). The units of the axes areall scaled by k_(max). There are two views so that the top and bottomsurfaces are both visible. The top and bottom surfaces are spherical(with different radii and centers), and the side surface is a rightcircular cone with its vertex at the origin.

In the limit of small bandwidth and low numerical aperture, the bandvolume shape approaches that of a circular cylinder. In this limit, theresolution in the axial direction is determined solely by the bandwidth,and the transverse resolution is determined by the numerical aperture,as is normally assumed in OCT. However, the band volume becomes lesscylindrical and more cone-shaped as the numerical aperture and bandwidthincrease, and axial and transverse resolutions are dependent on both thebandwidth and numerical aperture.

Simulation of Full-Field ISAM

A simulation is now presented to demonstrate inverse scattering infull-field OCT. For purposes of the simulation, an object (element 8shown in FIG. 3( a)) is taken as comprising a set of randomly placedpoint scatterers. The simulated object was imaged with a simulatedfull-field OCT system, and then the structure of the object wasreconstructed from the data. The simulated object volume cross-sectionalarea was 25 wavelengths in depth, and 40 by 40 wavelengths in thetransverse direction. The illumination source had a Gaussian spectrumwith a 40% fractional full-width-half-max bandwidth (corresponding, forexample, to 320 nm of bandwidth centered at 800 nm). The simulatednumerical aperture of the imaging objective was 0.5.

To synthesize the data, first the scattered field E_(s)(r′;k) wascalculated using Eq. (2), where the object η(r) was a collection ofrandomly chosen discrete points. From this, a two-dimensional Fouriertransform computed {tilde over (E)}_(s)(q;k). Then the synthesized datafunction was calculated by {tilde over (D)}(q;k)=A(k){tilde over(E)}_(s)(q;k){tilde over (P)}(q;k). Finally, a two-dimensional inverseFourier transform yielded D(r′;k). Eq. (10) was deliberately not used tocompute the data because using an alternate and more direct method ofcomputing the data provided a better test of the inverse scatteringmethod.

FIG. 3( c) shows the projection of the simulated data collapsed (summed)along one transverse direction. The units are in terms of the centerwavelength. Instead of showing the Fourier-transformed function {tildeover (D)}(r;k) itself, which would be difficult to interpret if it wasplotted directly, we show the inverse Fourier transform of {tilde over(D)}(r;k) with respect to k. It is the data on the focal plane arraythat would be observed if the delay τ were scanned, rather than thefrequency k, which is given by the intensity function of Eq. (12). Thefocus is on the top surface at zero depth, which also corresponds tozero delay. As can be seen, as the delay is increased from zero, thediffracted images of the point scatterers become increasingly large.This corresponds to the standard degradation in resolution one expectsfrom defocus when inverse scattering is not used.

To compute the image estimate η⁺(r) from the synthetic data D(r; k),first D(q;k) was computed using the two-dimensional Fourier transform.Next, Eq. (18) was used to compute {tilde over (η)}⁺{q+{circumflex over(z)}[k+k_(z)(q)]}. However, in practice to find η⁺(r) from {tilde over(η)}⁺(Q) numerically, one would likely use the three-dimensional inversediscrete Fourier transform. Unfortunately, Eq. (18) does not specify{tilde over (η)}⁺ in a form to which the inverse discrete Fouriertransform can be readily applied, because it is a function of the morecomplicated argument q+{circumflex over (z)}[k+k_(z)(q)]. In practice,this means that the discrete sampling of the function {tilde over (η)}⁺is uniform in the variables q and k and not in Q to which the inverseFourier transform can be directly applied. By using an interpolator, onecan compute the samples of {tilde over (η)}⁺ on points that are uniformin Q from existing samples of {tilde over (η)}^(†){q+{circumflex over(z)}[k+k_(z)(q)]} In this simulation, a one-dimensional cubic B-splineinterpolator was used to interpolate from the coordinates q+{circumflexover (z)}[k+k_(z)(q)] to Q. Because only the {circumflex over (z)}coordinate is not sampled uniformly, the resampling only needs to occuralong this direction.

Finally, after taking the three-dimensional inverse Fourier transform of{tilde over (η)}⁺(Q), the reconstruction η⁺(r) results, which is shownin FIG. 3( d). As can be seen, the reconstruction corrects for thediffraction of the data, and produces point-like images. FIG. 4 showsthree en face planes corresponding to the depths A, B, and C marked inFIG. 3( c). The left column is the time-domain data measured in each ofthe en face planes, and the right column is the image of the scattererscomputed by inverse scattering. Planes that are further from the focushave more diffuse images when viewed in the raw data because of theeffect of defocus. One can also see the interference between the imagesof adjacent scatterers. Despite the interference between scatterers,each point is clearly resolved with space-invariant resolution in thereconstructed image. This shows the algorithm correctly separates theinterference patterns from scatterers to produce high resolution images.

To show the overall improvement to the data, FIG. 5( a) is a volumeisosurface plot of the raw data, while the reconstructed computed imageis shown in FIG. 5( b). Again, the blurring of the data is increasinglyapparent with increasing distance from the focus plane at the top of thevolume. In addition, stripe-like features can be seen for theisosurfaces corresponding to interfering scatterers. This method cancorrect for the diffraction effects and produce point-like images inFIG. 5( b) for each of the scatterers. The planes of the scatterers neednot be so widely separated for the algorithm to distinguish them, butwas deliberately done to make the diffraction effects easier tovisualize.

There is an important difference in the reconstructions of full-fieldOCT and conventional scanned beam OCT. In conventional scanned beam OCT,it has been shown by Ralston et al., Inverse Scattering for OpticalCoherence Tomography, J. Opt. Soc. Am. A, vol. 23, pp. 1027-1037,(2006), incorporated herein by reference, that the magnitude of thesignal captured from scatterers away from the focus is inverselyproportional to the distance from the focus. In practice this places alimit on the axial range of the sample that can be imaged before thesignal-to-noise ratio becomes unacceptable. However, there is no suchattenuation of the signal away from the focus in the full-field OCTcase. The practical limit to the depth of full-field OCT is determinedby vignetting of the relayed field inside the relay telescope, and thescattering of the sample. However, this advantage may be offset becausefull-field OCT may be less able to discriminate between singlescattering and multiply scattered photons due to its multimodedetection.

Transverse-Scanned Focused Gaussian Beam Implementation

Preferred embodiments of the present invention implement computedinterferometric imaging by employing a fiber-optic Michelsoninterferometer 600 seeded by a source of femtosecond pulses, as nowdescribed with reference to FIG. 6. A spectral interferometer,designated generally by numeral 602, measures the interferometriccross-correlation between a fixed-delay reference pulse and a pulsereflected back from a sample 8. The measured spectrum on a line camera604 corresponds to the Fourier transform of the cross-correlationsignal, from which the amplitude and phase of the reflected field fromthe sample are inferred. The sample probe beam 606 is focused by amicroscope objective 608 so that the beam focus is at a fixed distanceinside the sample. At each position of the beam, the spectralinterferogram of the backscattered optical field is measured. The beamis laterally translated in two-dimensions through the sample by moving astage 610 or by steering the beam with a pair of galvanometer-drivenmirrors 612 before entering the objective.

Measurements are made using a femtosecond spectral interferometer 600.Fourier-domain, or frequency-swept, OCT, has been described by Choma etal., Sensitivity advantage of swept source and Fourier domain opticalcoherence tomography, Opt. Exr., vol. 111, pp. 2183-89 (2003), which isincorporated herein by reference. A femtosecond laser 614 (such assupplied by Kapteyn-Murnane Laboratories of Boulder, Colo.) deliversultrashort pulses to provide broadband illumination for the system. Inone embodiment of the invention, the center wavelength of the source is800 nm, with a bandwidth of 100 nm. These first-order field quantitiesfluctuate too rapidly to be detected directly, thus an opticalfiber-based Michelson interferometer is incorporated. The illuminationis divided by a 50:50 fiber-optic coupler (splitter) 616 between areference arm 618 containing a delay mirror 620 and a sample arm thatcontains a lens (objective) 622 to focus the Gaussian beam into thesample 8. Light returns from the sample and reference arms and isdirected into spectrometer 602. In the spectrometer, the light iscollimated with an achromatic lens 624 and dispersed off of a blazedgold diffraction grating 626, which, in one embodiment, has 830.3grooves per millimeter and a blaze angle of 19.70 degrees for a blazewavelength of 828 nm. To reduce lens aberrations, the dispersed opticalspectrum is focused using a pair of achromatic lenses 628. The focusedlight is incident upon a line-scan camera (L104k-2k, Basler, Inc.) 604which contains a 2048-element charge-coupled device (CCD) array with10×10 pm detection elements. Camera 604 operates at a maximum readoutrate of over 29 kHz, thus one axial scan can be captured during anexposure interval of about 34 μs. The data is sent to processor 630which may also govern operation of a galvanometric controller 642 andreceive a trigger derived therefrom in order to activate frameacquisition, for example.

We obtain a mathematical model of interferometric synthetic aperturemicroscopy (ISAM) by considering the propagation of the focused beamfrom the objective into the sample (into some volume V), scatteringwithin the sample (in the first Born approximation), the propagation ofthe scattered light back into the objective (over some surface Σ), andthe measurement of the cross-correlation with the reference pulse. Theexpression that models these steps (17) is given byS(r _(o) ,k)=A(k)∫_(Σ) d ² r∫ _(V) d ³ r′G(r′,r,k)g(r′−r _(o),k)η(r′)g(r−r _(o) ,k)  (21)where k is the wave number of the illumination, r₀ is the transverseposition of the Gaussian beam, g describes the normalized Gaussian beamprofile, A²(k) is the power spectral density of the source, G is theGreen function, and η is the susceptibility of the sample. Thenormalized beam profile is given by g(r,k)=W⁻²(k)e^(−r) ² ^(/2W) ²^((k))/2π, where W(k)=α/k, α=π/NA, and NA is the numerical aperture ofthe beam. The Green function is given by G(r′,r,k)=e^(ik(r−r′))/|r−r′|.After two-dimensional (2-D) Fourier transformation with respect to r₀,and further manipulation, the 2-D Fourier transform of the signal isgiven by the expression

$\begin{matrix}{{{\overset{\sim}{S}\left( {Q,k} \right)} = {{A(k)}{\int{{\mathbb{d}^{2}q}{\int{{\mathbb{d}z^{\prime}}\frac{{\mathbb{i}}\; 2\;\pi}{k_{z}(q)}{\mathbb{e}}^{{\mathbb{i}}\;{k_{z}{(q)}}{({z^{\prime} - z_{0}})}}{{\overset{\sim}{g}}_{0}\left( {q,k} \right)}{\mathbb{e}}^{{\mathbb{i}}\;{k_{z}{({Q - q})}}{({z^{\prime} - z_{0}})}}{{\overset{\sim}{g}}_{0}\left( {{Q - q},k} \right)}{\overset{\sim}{\eta}\left( {Q;z^{\prime}} \right)}}}}}}},} & (22)\end{matrix}$where k_(z)(q)=√{square root over (k²−a²)}, z₀ is the position of thebeam focus, {tilde over (η)}(Q,z′) is the 2-D transverse Fouriertransform of the susceptibility, and {tilde over (g)}₀(q,k)=e^(−q) ²^(α) ² ^(/2k) ² is the 2-D Fourier transform for the beam profile g(r,k)in the waist plane of the beam. After the expression for {tilde over(g)}₀ is substituted into Eq. (22), and an asymptotic expansion of{tilde over (S)} is made for large α², this relationship reduces to

$\begin{matrix}{{{\overset{\sim}{S}\left( {Q,k} \right)} = {{A(k)}\left( {\frac{{\mathbb{i}}\; 2\;\pi^{2}}{k_{z}\left( {Q/2} \right)}\frac{k^{2}}{\alpha^{2}}{\mathbb{e}}^{{- 2}\;{\mathbb{i}}\;{k_{z}{({Q/2})}}z_{0}}{\mathbb{e}}^{- \frac{\alpha^{2}Q^{2}}{4k^{2}}}} \right){\overset{\sim}{\eta}\left( {Q;{{- 2}{k_{z}\left( {Q/2} \right)}}} \right)}}},} & (23)\end{matrix}$where {tilde over ({tilde over (η)} is the 3-D Fourier transform of η,i.e. the one dimensional Fourier transform of {tilde over (η)}(Q;z) withrespect to z. This expansion is valid even when NA≈1 because α² issufficiently large for the first term of the expansion to dominate. Eq.(23) relates the 3-D Fourier transform of the object susceptibility tothe 2-D Fourier transform of the signal. Implicit in this formula is adiagonal linear integral operator in the 3-D Fourier space of thesusceptibility, and so the achievable resolution is spatially-invariantand does not depend on the proximity to the focus.

Because of the simple relationship between the susceptibility and thesignal, ISAM can be implemented efficiently by resampling orinterpolating the data in a manner analogous to the numericalimplementation of the Fourier projection slice theorem, as described inNatterer, The Radon Transform, (Wiley, 1986, incorporated herein byreference) and as used in x-ray computed tomography or syntheticaperture radar (SAR), but the resampling grid for ISAM is hyperbolicrather than polar. In addition, since Eq. (23) is a multiplicative (ordiagonal) form, generalization to a regularized inverse method such asTikhonov regularization (Tikhonov, Dokl. Akad. Nauk SSR, vol. 39, p.195, 1943) is straightforward.

Regardless of the type of detection used, the beam orientation is fixedwhile the axis is translated over the sample on a Cartesian grid indirections orthogonal to the beam axis and subsequent axial scans aredisplayed on adjacent lines to form a tomogram. Suppose the beam isgenerated by light emerging from an optical fiber and then propagatingthrough a lens to form a focused beam incident on the sample. With theaxis of the fiber passing through the point r₀ in the z=0 plane and withthe waist plane of the focused field at z=z₀, the incident field may bedescribed by the power spectrum |A(k)|² and a normalized mode g suchthatU ₁(r,r ₀ ,k)=A(k)g(r−r ₀).  (24)The beam may be described in a plane wave decomposition,

$\begin{matrix}{{{g\left( {{r - r_{0}},k} \right)} = {\frac{1}{\left( {2\;\pi} \right)^{2}}{\int{{\mathbb{d}^{2}q}\;{\mathbb{e}}^{{\mathbb{i}}\;{q \cdot {({r - r_{0}})}}}{\mathbb{e}}^{{\mathbb{i}}\;{k_{z}{(q)}}{({z - z_{0}})}}{{\overset{\sim}{g}}_{0}\left( {q,k} \right)}}}}},} & (25)\end{matrix}$

where {tilde over (g)}₀ is the two dimensional Fourier transform of g inthe z=z₀ plane, and the dispersion relation is given by k_(z)=√{squareroot over (k²−q²)}. The beam waist is assumed to depend on the wavenumber, as W₀(k)=α/k where α=π/NA and NA is the numerical aperture ofthe output lens. Thus{tilde over (g)} ₀(q,k)=e ^(−q) ² ^(W) ⁰ ² ^(/2) =e ^(−q) ² ^(α) ²^(/(2k) ³ ⁾.  (26)

The scattered field, within the first Born approximation, is given byU _(s)(r,r ₀ ,k)=∫d ³ r′G(r′,r,k)U ₁(r′,r ₀ ,k)η(r′).  (27)Making use of the expressions above for the incident field,U _(s)(r,r ₀ ,k)=A(k)∫d ³ r′G(r′,r,k)g(r′−r ₀ ,k)η(r′).  (28)The signal coupled back in to the fiber is given by the projection ofthe backscattered field onto the fiber mode g at the exit plane of thefiber. ThusS(r ₀ ,k)=∫_(z=0) d ² rU(r,r ₀ ,k)g(r−r ₀ ,k),  (29)which becomesS(r ₀ ,k)=A(k)∫_(z=0) =d ² r∫d ³ r′G(r′,r,k)g(r′−r ₀ ,k)g(r−r ₀,k)η(r′).  (30)The Green function for free-space is given by the angular spectrum

$\begin{matrix}{{{G\left( {r^{\prime},r,k} \right)} = {\frac{\mathbb{i}}{2\;\pi}{\int{{\mathbb{d}^{2}q}\;{\mathbb{e}}^{{\mathbb{i}}\;{q \cdot {({r - r^{\prime}})}}}\frac{{\mathbb{e}}^{{- {\mathbb{i}}}\;{k_{z}{(q)}}{({z - z^{\prime}})}}}{k_{z}(q)}}}}},} & (31)\end{matrix}$where it is assumed that the scatterers are all located such that z<z′for the whole support of the scatterers. Making use of this expressionand Eq. (30), it may be seen that the two-dimensional Fourier transformof S with respect to r₀ is given by the expression

$\begin{matrix}{{\overset{\sim}{S}\left( {Q,k} \right)} = {{\mathbb{i}}\; 2\;\pi\;{A(k)}{\int{{\mathbb{d}^{2}q}{\int{{\mathbb{d}z^{\prime}}\frac{1}{k_{z}(q)}{\mathbb{e}}^{{\mathbb{i}}\;{k_{z}{(q)}}{({z^{\prime} - z_{0}})}}{\mathbb{e}}^{{\mathbb{i}}\;{k_{z}{({q - Q})}}{({z^{\prime} - z_{0}})}}{\mathbb{e}}^{\frac{{- \alpha^{2}}q^{2}}{2k^{2}}}{\mathbb{e}}^{\frac{{- \alpha^{2}}{{q - Q}}^{2}}{2k^{2}}}{{\overset{\sim}{\eta}\left( {Q,z^{\prime}} \right)}.}}}}}}} & (32)\end{matrix}$This equation may be solved for η by blunt numerical methods. Suchmethods are numerically expensive. An analytic result may be obtained byconsidering the shifted form of the integral

$\begin{matrix}{{\overset{\sim}{S}\left( {Q,k} \right)} = {{\mathbb{i}}\; 2\;\pi\;{A(k)}{\int{{\mathbb{d}^{2}q}{\int{{\mathbb{d}z^{\prime}}\frac{1}{k_{z}(q)}{\mathbb{e}}^{{\mathbb{i}}\;{k_{z}{(q)}}{({z^{\prime} - z_{0}})}}{\mathbb{e}}^{{\mathbb{i}}\;{k_{z}{({q - Q})}}{({z^{\prime} - z_{0}})}}{\mathbb{e}}^{\frac{{- \alpha^{2}}Q^{2}}{4k^{2}}}{\mathbb{e}}^{\frac{{- \alpha^{2}}{{q - {Q/2}}}^{2}}{k^{2}}}{{\overset{\sim}{\eta}\left( {Q,z^{\prime}} \right)}.}}}}}}} & (33)\end{matrix}$For large values of α this integral may be evaluated asymptotically. Theintegrand, modulo the Gaussian, may be expanded in a Taylor seriesaround the point q=Q/2,

$\begin{matrix}{\frac{{\mathbb{e}}^{{{\mathbb{i}}{\lbrack{{k_{z}{(q)}} + {k_{z}{({q - {Q/2}})}}}\rbrack}}{({z_{0} - z^{\prime}})}}}{k_{z}(q)} = \left. {\frac{{\mathbb{e}}^{2\;{\mathbb{i}}\;{k_{z}{({Q/2})}}{({z^{\prime} - z_{0}})}}}{k_{z}\left( {Q/2} \right)} + {q \cdot {\nabla_{q}\frac{{\mathbb{e}}^{{{\mathbb{i}}{\lbrack{{k_{z}{(q)}} + {k_{z}{({q - {Q/2}})}}}\rbrack}}{({z_{0} - z^{\prime}})}}}{k_{z}(q)}}}} \middle| {}_{q = {Q/2}}{+ \ldots} \right.} & (34)\end{matrix}$Replacing this part of the integrand, the leading term is given by anintegral over the constant term in the Taylor expansion:

$\begin{matrix}{{\overset{\sim}{S}\left( {Q,k} \right)} = {{\mathbb{i}}\; 2\;\pi\;{A(k)}{\mathbb{e}}^{\frac{{- \alpha^{2}}Q^{2}}{4k^{2}}}{\int{{\mathbb{d}z^{\prime}}\frac{{\mathbb{e}}^{2\;{\mathbb{i}}\;{k_{z}{({Q/2})}}{({z^{\prime} - z_{0}})}}}{k_{z}\left( {Q/2} \right)}{\int{{\mathbb{d}^{2}q}\;{\mathbb{e}}^{\frac{{- \alpha^{2}}{{q - {Q/2}}}^{2}}{k^{2}}}{{\overset{\sim}{\eta}\left( {Q,z^{\prime}} \right)}.}}}}}}} & (35)\end{matrix}$The Gaussian integral may be easily carried out and the remainingintegral is seen to be a Fourier transform with respect to z′,

$\begin{matrix}{{{\overset{\sim}{S}\left( {Q,k} \right)} = {\frac{k^{2}}{\alpha^{2}}{\mathbb{i}}\; 2\pi^{2}{A(k)}\frac{{\mathbb{e}}^{{- 2}{\mathbb{i}}\;{k_{z}{({Q/2})}}z_{0}}}{k_{z}\left( {Q/2} \right)}{\mathbb{e}}^{\frac{{- \alpha^{2}}Q^{2}}{4\; k^{2}}}{\overset{\overset{\sim}{\sim}}{\eta}\left\lbrack {Q,{{- 2}{k_{z}\left( {Q/2} \right)}}} \right\rbrack}}},} & (36)\end{matrix}$where {tilde over ({tilde over (η)} is the three-dimensional Fouriertransform of η. The next term in the expansion yields a contributionproportional to α⁻⁴. In the extreme limit that NA→1, it may be seen thatα→π so that we expect the leading term approximation to be sufficienteven in the case of high numerical aperture. It might be noted that thisexpansion is distinct from the paraxial approximation (essentially asmall |q| expansion of k_(z)(q)) which fails as NA→1. Eq. (36) expressesa resampling scheme illustrated in FIG. 1. To find an appropriateregularization scheme, we will write

$\begin{matrix}{{{\overset{\sim}{S}\left( {Q,k} \right)} = {\int{{\mathbb{d}\beta}\;{H\left( {Q,k,\beta} \right)}{\overset{\overset{\sim}{\sim}}{\eta}\left( {Q,\beta} \right)}}}},{where}} & (37) \\\begin{matrix}{{H\left( {Q,k,\beta} \right)} = {\frac{k^{2}}{\alpha^{2}}{\mathbb{i}}\; 2\pi^{2}{A(k)}\frac{{\mathbb{e}}^{{- 2}{\mathbb{i}}\;{k_{z}{({Q/2})}}z_{0}}}{k_{z}\left( {Q/2} \right)}{\mathbb{e}}^{\frac{{- \alpha^{2}}Q^{2}}{4\; k^{2}}}{\delta\left\lbrack {\beta + {2\;{k_{z}\left( {Q/2} \right)}}} \right\rbrack}}} \\{\equiv {{f\left( {Q,k,\beta} \right)}{{\delta\left\lbrack {\beta + {2\;{k_{z}\left( {Q/2} \right)}}} \right\rbrack}.}}}\end{matrix} & (38)\end{matrix}$Then the kernel of the normal operator is given by the expression

$\begin{matrix}{{H^{*}{H\left( {Q,\beta,\beta^{\prime}} \right)}} \equiv {{{f\left( {Q,{{1/2}\sqrt{\beta^{2} + Q^{2}}},\beta} \right)}}^{2}\frac{\beta}{2\sqrt{\beta^{2} + Q^{2}}}{{\delta\left( {\beta - \beta^{\prime}} \right)}.}}} & (39)\end{matrix}$And the kernel of the Tikhonov regularized psuedo-inverse, with whitenoise N is given by the expression

$\begin{matrix}{{H^{+}\left( {Q,{k;\beta}} \right)} = {\frac{{f^{*}\left( {Q,k,\beta} \right)}{\delta\left( {k - {{1/2}\sqrt{\beta^{2} + Q^{2}}}} \right)}}{{{f\left( {Q,k,\beta} \right)}}^{2} + {2{{Nk}/{k_{z}\left( {Q/2} \right)}}}}.}} & (40)\end{matrix}$The object structure is then given by

$\begin{matrix}{{{\overset{\overset{\sim}{\sim}}{\eta}}^{+}\left( {Q,\beta} \right)} = \left\lbrack \frac{{f^{*}\left( {Q,k,\beta} \right)}{\overset{\sim}{S}\left( {Q,k} \right)}}{{{f\left( {Q,k,\beta} \right)}}^{2} + {2\;{{Nk}/{k_{z}\left( {Q/2} \right)}}}} \right\rbrack_{k = {\frac{1}{2}\sqrt{\beta^{2} + Q^{2}}}}} & (41)\end{matrix}$The object structure in the coordinate domain is obtained by applyingthe three-dimensional inverse Fourier transform.

To achieve phase stability of the signal, a microscope coverslip (notshown) may be placed on top of sample 8 and the top reflection from theair-coverslip interface acts as a fixed reference delay relative to theobject. The delay fluctuations of the interferometer are removed fromeach cross-correlation interferogram by locating the air-coverslipreflection in each interferogram, estimating the phase and group delayof the reflection, and applying the opposite phase and group delay tothe entire interferogram. Prior to processing, the spectra, eachrepresenting a column of depth-dependent data, are assembled adjacentlyas the beam is transversely scanned over the sample. The detecteddigital signal is interpolated to account for the non-uniform samplingof the spectrum and to compensate for up to third-order dispersion.Specifically, the signal is interpolated by a factor of two by aband-limiting interpolator implemented using the fast Fourier transform(FFT). This prepares the signal for the cubic B-spline interpolation,which has a transfer function with an amplitude that attenuatesfrequencies close to the Nyquist limit. The cubic B-spline interpolatorresamples the spectrum to a uniform frequency space according to acalibration procedure utilizing a single reflector placed at the focusof the objective lens. Sample movement, inducing phase and group delaychanges, is tracked using a reference microscope coverslip, and thedeviations are corrected. At this point, the quantity S(r₀, k), in Eq.(1) has been estimated. Next, the two-dimensional FFT in the transversedirections is calculated to yield {tilde over (S)}(Q,k). Then, thenon-uniform ISAM resampling and filtering of Eq. (3) using cubicB-splines is implemented to yield {tilde over ({tilde over (η)}.Finally, the 3-D inverse FFT is used to attain the ISAM reconstruction,an estimate of η(r).

By application of ISAM techniques as described, in vivo imaging mayadvantageously be performed on larger volumes of tissue than volumesthat would otherwise have to be resected. Furthermore, ISAM achieveshigh-speed, high-resolution imagery without need for the timelyprocessing, sectioning, and staining of a resection.

With the use of near-infrared light, high-resolution ISAM facilitatesthe noninvasive monitoring of cellular and nuclear development withpenetration depths up to 3 mm. Of course, in regions of extremescattering or absorption penetration depths may be reduced.

Image formation algorithms are characterized differently than imagepost-processing routines. In particular, ISAM is a unique imageformation method that utilizes the fundamental resolution capabilitiesof the acquired optical signal based on the physics of the scatteringwithin the detection beam. In contrast, a simple image post-processingalgorithm may attempt to extrapolate beyond the information inherent animage. For example, maximum likelihood, blind-deconvolution, orentropy-constrained algorithms can effectively produce energy compactionas an estimated solution to an image reconstruction. Such estimation mayincorrectly reconstruct features in an image, thus misrepresenting thetrue biological structure and potentially leading to an incorrectdiagnosis. ISAM is not based on estimation, thus such misrepresentationsdo not exist. Furthermore, estimation algorithms often exceed thelimited computational complexity necessary for real-time imaging. TheISAM image formation algorithm can be implemented with computationalcomplexity of O(NlogN), where N is the number of volume elements toresolve, which makes ISAM amenable to real-time imaging. Furthermore,the ISAM algorithm can be applied to planes as well as volumes, thusenhancing cross-sectional imaging.

Azimuthally-Scanned Implementation

Embodiments of the invention are now described in which a focused beamis directed perpendicular to an OCT catheter, which might, for example,be incorporated into an endoscope. An endoscope, used for exploring thetubular lumens within the human gastrointestinal tract, typicallyconsists of a long, flexible device of 1 cm diameter or less. Inside theendoscope, in addition to a source of white light illumination andoptics for imaging, for example, on a charge-coupled device (CCD)detector, are working channels through which various instruments forbiopsy or manipulating tissue are passed. For example, tissue biopsysamples can be extracted and withdrawn by forceps or suction. Smallerdiameter catheters are used in the cardiovascular system, e.g. for theinsertion of balloons for angioplasty or to deploy stents. Intravascularcatheters minimize invasiveness and provide access to vascular lesionsassociated with cardiovascular disease.

An azimuthally-scanned OCT system would typically include a workingchannel containing a single-mode optical fiber, a focusing lens(typically a graded index lens cemented or fused to the fiber), and aright-angle prism or a custom cleaved surface for reflecting the beam by90 degrees to the side of the catheter.

Fiber-optic OCT catheters have been integrated with endoscopes to imagethe esophagus, colon, and other internal organs and mucosal tissue, asdescribed, for example, by Tearney et al., In vivo endoscopic opticalbiopsy with optical coherence tomography, Science, vol. 276, pp. 2037-39(1997), incorporated herein by reference. In instruments based onfiber-optic OCT, the illumination originates inside the object ortubular lumen being imaged, and is usually scanned azimuthally aroundthe long axis of the catheter. As the catheter is azimuthally scannedand translated along the long-axis of the catheter, a 3-D image of theobject is acquired. Because the beam is typically focused at a fixeddistance from the catheter, the depth-of-focus of the resulting imagesis confined to a narrow annulus.

By rotating the catheter about its long-axis, the beam may be directedalong any path perpendicular to the axis. By pushing or pulling thecatheter, the beam is translated along the long-axis of the catheter.Together these two degrees of freedom enable the instrument to scan acylindrically shaped volume around the catheter. Typical imaging withthis catheter design involves acquisition of axial scans (either in thetime or frequency domain) while rotating the catheter through 360degrees, advancing the catheter a small distance along its long-axis,and repeating the measurement. After acquisition, one possesses a dataset parameterized by the illumination frequency (or time delay), theangular coordinate of the catheter during the scan, and thetranslational position of the catheter along its axis. With our solutionof the inverse problem, we infer the object susceptibility from thesedata.

An algorithm that infers the susceptibility of a scatterer from thesignal acquired in angularly scanned OCT is now described. These may beadvantageously employed in catheter-based optical coherence tomography,but the scope of the present invention is not limited and may includeapplication in endoscopic or intravascular ultrasound as well. Otherapplications may include acoustic, sonar, and seismic sensing where theimaged object is close to a focused transducer, and radar sensing ofobjects near a rotating dish antenna.

The Forward Problem for Azimuthally-Scanned ISAM

In contradistinction to the forgoing discussion wherein the illuminatingGaussian beam was translated in a straight line perpendicular to itsaxis, in the following discussion, rotation of the Gaussian beam isconsidered about the origin.

We consider an experiment in which a Gaussian beam originates at aposition with Cartesian coordinates (0, p, 0). Let us denote Cartesiancoordinates fixed relative to the sample by r=(x, y, z) and let usdenote Cartesian coordinates fixed relative to the beam by r′=(x′, y′,z′). For each fixed axial position of the fiber y=y′=p. The beam isdirected at an angle θ from the z-axis, and along the z′ axis. Thecoordinates may be related by a rotation matrix R(θ) so that r=R(θ)r′where

$\begin{matrix}{{R(\theta)} = {\begin{pmatrix}{\cos\;\theta} & 0 & {\sin\;\theta} \\0 & 1 & 0 \\{{- \sin}\;\theta} & 0 & {\cos\;\theta}\end{pmatrix}.}} & (42)\end{matrix}$

The beam is focused a distance z₀ from they axis. The field ispolychromatic with power spectrum A²(k) where k=ω/c is the wave numberassociated with frequency ω. The width of the beam waist is a functionof frequency given by W(k)=α/k, where α=π/NA, and NA, as above, is thenumerical aperture of the focused beam. The beam is rotationally scannedso that the signal is sampled for all angles −π≦θ<π, and the beam isalso translated along they axis to cover all axial coordinates p. FIG. 7illustrates this notation.

In the discussion above, it was assumed that the direction ofpropagation of the beam was fixed to be along the z direction. Thelocation of the center of the beam in the waist plane was given by r₀.Then the signal measured in the interferometer is given by theexpression {tilde over (S)}(r,k), which is given by

$\begin{matrix}{{{\overset{\sim}{S}\left( {r,k} \right)} = {{{\mathbb{i}}\left( {2\pi} \right)}^{- 2}{A(k)}k^{- 1}{\int_{V}^{\;}\ {{\mathbb{d}^{3}r}\;\eta\;(r){f^{2}\left( {{r - r_{0}};k} \right)}}}}},} & (43)\end{matrix}$Where η(r) is the susceptibility of the sample being probed, f² (r′;k)is given by the expression:

$\begin{matrix}\begin{matrix}{{f^{2}\left( {r^{\prime};k} \right)} = {\frac{1}{\left( {2\;\pi} \right)^{2}}{\int{{\mathbb{d}^{2}{{\xi exp}\left( {{- {\mathbb{i}\xi}} \cdot r^{\prime}} \right)}}\frac{1}{2}\left( {\frac{\alpha^{2}}{k^{2}} + \frac{{\mathbb{i}}\left( {z^{\prime} - z_{0}} \right)}{k}} \right)^{- 1}}}}} \\{{\exp\left( {- \frac{\xi^{2}\alpha^{2}}{4\; k^{2}}} \right)}{\exp\;\left\lbrack {{{\mathbb{i}}\left( {z^{\prime} - z_{0}} \right)}\sqrt{\left( {2\; k} \right)^{2} - \xi^{2}}} \right\rbrack}}\end{matrix} & (44)\end{matrix}$where ξ=(ξ_(x),ξ_(y),0) and the integral is over the ξ_(x), ξ_(y) plane.Note that we do not now make the paraxial approximation for the phaseterm. The signal depends on frequency, position along the y-axis, andthe angle of propagation of the beam as described above. This signal,{tilde over (S)}(k,p,θ), may be found from Eq. (39) by writing theintegrand in the coordinates stationary with respect to the beam. Thuswe obtain

$\begin{matrix}{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {{{\mathbb{i}}\left( {2\;\pi} \right)}^{- 2}A(k)k^{- 1}{\int_{V}^{\;}\ {{\mathbb{d}^{3}r^{\prime}}{\eta\left\lbrack {{R(\theta)}r^{\prime}} \right\rbrack}{{f^{2}\left( {{r^{\prime} - {p\hat{y}}};k} \right)}.}}}}} & (45)\end{matrix}$

By substituting Eq. (44) into Eq. (45) and rearranging terms, we find

$\begin{matrix}\begin{matrix}{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {\frac{\mathbb{i}}{2}{A(k)}k^{- 1}{\int{{\mathbb{d}^{2}{{\xi exp}\left\lbrack {{- {\mathbb{i}}}\; z_{0}\sqrt{\left( {2k} \right)^{2} - \xi^{2}}} \right\rbrack}}{\exp\left( {- \frac{\xi^{2}\alpha^{2}}{4\; k^{2}}} \right)}}}}} \\{\int{{\mathbb{d}^{3}r^{\prime}}{\exp\left\lbrack {{- {\mathbb{i}}}\;{\xi \cdot \left( {r^{\prime} - {p\hat{y}}} \right)}} \right\rbrack}{\eta\left\lbrack {{R(\theta)}r^{\prime}} \right\rbrack}}} \\{\left\lbrack {\frac{\alpha^{2}}{k^{2}} + \frac{{\mathbb{i}}\left( {z^{\prime} - z_{0}} \right)}{k}} \right\rbrack^{- 1}{{\exp\left\lbrack {{\mathbb{i}}\; z^{\prime}\sqrt{\left( {2\; k} \right)^{2} - \xi^{2}}} \right\rbrack}.}}\end{matrix} & (46)\end{matrix}$

In our analysis of the OCT inverse problem on a Cartesian grid, we foundthat under certain reasonable approximations, the data could be relatedto the object susceptibility through a resampling scheme in the Fourierdomain. We derive a similar relation here. To do so, it will beadvantageous to replace

$\left\lbrack {\frac{\alpha^{2}}{k^{2}} + \frac{{\mathbb{i}}\left( {z^{\prime} - z_{0}} \right)}{k}} \right\rbrack$with an approximation commensurate with the natural geometry of theproblem. Explicitly, we replace z′ with ρ′=√{square root over(z′²+x′²)}. For most OCT systems, the bandwidth is a small fraction ofthe central frequency and so we replace

$\frac{1}{k^{2}}{with}{\frac{1}{k\; k_{0}}.}$Thus the factor

$\left\lbrack {\frac{\alpha^{2}}{k^{2}} + \frac{{\mathbb{i}}\left( {z^{\prime} - z_{0}} \right)}{k}} \right\rbrack$is replaced by

${\frac{1}{k}\left\lbrack {\frac{\alpha^{2}}{k_{0}} + {{\mathbb{i}}\left( {\sqrt{x^{\prime^{2}} + z^{\prime^{2}}} - z_{0}} \right)}} \right\rbrack}.$This expression is slowly varying relative to the rapidly varying phaseof the term exp└iz′√{square root over ((2k)²−ξ²)}┘, and soapproximations to it tend not to change the result greatly. With thissubstitution,

$\begin{matrix}\begin{matrix}{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {\frac{\mathbb{i}}{2}{A(k)}{\int{{\mathbb{d}^{2}{{\xi exp}\left\lbrack {{- {\mathbb{i}}}\; z_{0}\sqrt{\left( {2k} \right)^{2} - \xi^{2}}} \right\rbrack}}{\exp\left( {- \frac{\xi^{2}\alpha^{2}}{4\; k^{2}}} \right)}}}}} \\{\int{{\mathbb{d}^{3}r^{\prime}}{\exp\left\lbrack {{- {\mathbb{i}}}\;{\xi \cdot \left( {r^{\prime} - {p\hat{y}}} \right)}} \right\rbrack}{{\eta\left\lbrack {{R(\theta)}r^{\prime}} \right\rbrack}\left\lbrack {\frac{\alpha^{2}}{k_{0}} + {{\mathbb{i}}\left( {\rho^{\prime} - z_{0}} \right)}} \right\rbrack}^{- 1}}} \\{{\exp\left\lbrack {{\mathbb{i}}\; z^{\prime}\sqrt{\left( {2\; k} \right)^{2} - \xi^{2}}} \right\rbrack}.}\end{matrix} & (47)\end{matrix}$

To evaluate this integral, we change variables in the inner integral tothe coordinates stationary in the reference frame of the sample,

$\begin{matrix}{{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {\frac{\mathbb{i}}{2}{A(k)}{\int{{\mathbb{d}^{2}{{\xi exp}\left\lbrack {{- {\mathbb{i}}}\; z_{0}\sqrt{\begin{matrix}{\left( {2k} \right)^{2} -} \\\xi^{2}\end{matrix}}} \right\rbrack}}{\exp\left( {- \frac{\xi^{2}\alpha^{2}}{4k^{2}}} \right)}{\exp\left( {{\mathbb{i}}\;\xi_{y}p} \right)}}}}}{{\int{{\mathbb{d}^{3}r}\;\exp\left\{ {{- {{\mathbb{i}}\left\lbrack {\xi - {\hat{z}\sqrt{\left( {2k} \right)^{2} - \xi^{2}}}} \right\rbrack}} \cdot {R\left\lbrack {\left( {- \theta} \right)r} \right\rbrack}} \right\}{{\eta(r)}\left\lbrack {\frac{\alpha^{2}}{k_{0}} + {{\mathbb{i}}\left( {\rho - z_{0}} \right)}} \right\rbrack}^{- 1}}},}} & (48)\end{matrix}$where ρ=ρ′=√{square root over (x²+z²)}. It may be seen that the integralover r results in a Fourier transform if we note that └ξ−{circumflexover (z)}√{square root over ((2k)²−ξ²)}┘·R(−θ)r=R(θ)└ξ−{circumflex over(z)}√{square root over ((2k)²−ξ²)}┘·r, after which we obtain

$\begin{matrix}{{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {\frac{\mathbb{i}}{2}{A(k)}{\int{{\mathbb{d}^{2}{{\xi exp}\left\lbrack {{\mathbb{i}}\; z_{0}\sqrt{\left( {2k} \right)^{2} - \xi^{2}}} \right\rbrack}}{\exp\left( {- \frac{\xi^{2}\alpha^{2}}{4k^{2}}} \right)}{\exp\left( {{\mathbb{i}}\;\xi_{y}p} \right)}}}}}{\overset{\sim}{\eta}\left\{ {- {{R(\theta)}\left\lbrack {\xi - {\hat{z}\sqrt{\left( {2k} \right)^{2} - \xi^{2}}}} \right\rbrack}} \right\}}} & (49)\end{matrix}$where {tilde over (η)}(β) is the weighted Fourier transform of η(r)given by

$\begin{matrix}{{\eta(\beta)} = {\int{{\mathbb{d}^{3}r}\;{\exp\left( {{\mathbb{i}}\;{r \cdot \beta}} \right)}{{{\eta(r)}\left\lbrack {\frac{\alpha^{2}}{k_{0}} + {{\mathbb{i}}\left( {\rho - z_{0}} \right)}} \right\rbrack}^{- 1}.}}}} & (50)\end{matrix}$

To change the integral over ξ to a cyclic convolution, we make thesubstitution √{square root over ((2k)²−ξ_(y) ²)} cos φ=ξ_(x) so that√{square root over ((2k)²−ξ_(y) ²)} sin φ=√{square root over((2k)²−ξ²)}, after which we obtain

$\begin{matrix}{{{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {\frac{\mathbb{i}}{2}{A(k)}{\int{{\mathbb{d}\xi_{y}}{\exp\left( {{\mathbb{i}}\;\xi_{y}p} \right)}}}}}\mspace{20mu}{\int_{0}^{\pi}{\mathbb{d}\left\{ {\left\lbrack {\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}\sin\;\phi} \right\rbrack{\exp\left\lbrack {{- {\mathbb{i}}}\; z_{0}\sqrt{\begin{matrix}{\left( {2k} \right)^{2} -} \\\xi^{2}\end{matrix}}} \right\rbrack}{\exp\left( {- \frac{\xi^{2}\alpha^{2}}{4k^{2}}} \right)}} \right\}}}{\overset{\sim}{\eta}{\left\{ {- {{R(\theta)}\left\lbrack {{\hat{x}\cos\;\phi\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}} + {\hat{y}\xi_{y}} - {\hat{z}\sin\;\phi\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}}} \right\rbrack}} \right\}.}}}\;} & (51)\end{matrix}$

For brevity, we define the kernel function K(k,ξ_(y),φ).

$\begin{matrix}{{K\left( {k,\xi_{y},\phi} \right)} = {\frac{\mathbb{i}}{2}{{A(k)}\left\lbrack {\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}\sin\;\phi} \right\rbrack}{\exp\left\lbrack {{- {\mathbb{i}}}\; z_{0}\sqrt{\left( {2k} \right)^{2} - \xi^{2}}} \right\rbrack}{{\exp\left( {- \frac{\xi^{2}\alpha^{2}}{4k^{2}}} \right)}.}}} & (52)\end{matrix}$

We note that the cos φ next to x and the sin φ next to {circumflex over(z)} in Eq. (51) effect a rotation in the x-z plane through an angle −φof a vector x√{square root over ((2k)²−ξ_(y) ²)}. Given this, we canexpress Eq. (51) as a cyclic convolution:

$\begin{matrix}{{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {\int{{\mathbb{d}\xi_{y}}{\exp\left( {{\mathbb{i}}\;\xi_{y}p} \right)}}}}{\int_{0}^{\pi}{{\mathbb{d}\phi}\;{K\left( {k,\xi_{y},\phi} \right)}\overset{\sim}{\eta}{\left\{ {- {{R\left( {\theta - \phi} \right)}\left\lbrack {{\hat{x}\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}} + {\hat{y}\xi_{y}}} \right\rbrack}} \right\}.}}}} & (53)\end{matrix}$

By combining the rotations R(θ) and R(−φ), we find the integral over φis a cyclic convolution. This cyclic convolution can be performedefficiently using the product of Fourier series. To put Eq. (51) intodiagonal form, we define the following functions of the data, thekernel, and the structure function:

$\begin{matrix}{{{\overset{\overset{\sim}{\sim}}{S}\left( {k,\xi_{p},n_{\theta}} \right)} = {\int_{- \infty}^{\infty}{\int_{- \pi}^{\pi}{{\mathbb{d}p}{\mathbb{d}{{\theta exp}\left( {{\mathbb{i}}\; p\;\xi_{p}} \right)}}{\exp\left( {{\mathbb{i}}\;\theta\; n_{\theta}} \right)}{\overset{\sim}{S}\left( {k,p,\theta} \right)}}}}},} & (54) \\{{{\overset{\sim}{K}\left( {k,\xi_{y},n_{\theta}} \right)} = {\int_{0}^{\pi}{{\mathbb{d}\theta}\;{\exp\left( {{\mathbb{i}}\;\theta\; n_{\theta}} \right)}{K\left( {k,\xi_{y},\theta} \right)}}}},} & (55) \\{{\overset{\overset{\sim}{\sim}}{\eta}\left( {k,\xi_{y},n_{\theta}} \right)} = {\int_{- \pi}^{\pi}{{\mathbb{d}\theta}\;{\exp\left( {{\mathbb{i}}\;\theta\; n_{\theta}} \right)}\overset{\sim}{\eta}{\left\{ {- \left\lbrack {{\hat{x}\cos\;\theta\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}} + {\hat{y}\xi_{y}} + {\hat{z}\sin\;\theta\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}}} \right\rbrack} \right\}.}}}} & (56)\end{matrix}$Where n_(θ) is an integer on [−∞, ∞]. If we insert the definitions ofEqs. (36)-(38) into Eq. (35), we find the following relationship:{tilde over ({tilde over (S)}(k,ξ _(p) ,n _(θ))={tilde over (K)}(k,ξ_(p) ,n _(θ)){tilde over ({tilde over (η)}(k,ξ _(p) ,n _(θ)).  (57)

In this form we see that {tilde over ({tilde over (S)} and {tilde over({tilde over (η)} are related by a diagonal integral operator whosekernel is {tilde over (K)}(k′,ξ_(p′),n_(θ′))δ(k−k′)δ(ξ_(p)−ξ_(p′))/δ_(n)_(θ) ,n_(θ′). Explicitly S=K{tilde over ({tilde over (η)} where{circumflex over (K)} is the integral operator

$\begin{matrix}{{\left\lfloor {K\;\overset{\overset{\sim}{\sim}}{\eta}} \right\rfloor\left( {k,\xi_{p},n_{\theta}} \right)} = {\int{{\mathbb{d}k^{\prime}}{\int{{\mathbb{d}\xi_{p^{\prime}}}{\sum\limits_{n_{\theta}}{{\overset{\sim}{K}\left( {k^{\prime},\xi_{p^{\prime}},n_{\theta^{\prime}}} \right)}{\delta\left( {k - k^{\prime}} \right)}{\delta\left( {\xi_{p} - \xi_{p^{\prime}}} \right)}\delta_{n_{\theta},n_{\theta^{\prime}}}{{\overset{\overset{\sim}{\sim}}{\eta}\left( {k^{\prime},\xi_{p^{\prime}},n_{\theta^{\prime}}} \right)}.}}}}}}}} & (58)\end{matrix}$

This diagonal operator will simplify finding solutions to specificinverse problems.

The Inverse Problem for Azimuthally-Scanned ISAM

Eq. (57) defines a linear relationship between the object structure anddata. To better understand how to invert this relationship, therelationship between the data {tilde over (S)}(k, p, θ) and the objectη(r) is written explicitly:

$\begin{matrix}{{{\overset{\sim}{S}\left( {k,p,\theta} \right)} = {\frac{1}{4\;\pi^{2}}{\int{{\mathbb{d}\xi_{y}}{\exp\left( {{- {\mathbb{i}}}\;\xi_{y}p} \right)}{\sum\limits_{n_{\theta} = {- \infty}}^{\infty}{{\exp\left( {{- {\mathbb{i}}}\;\theta\; n_{\theta}} \right)}{K\left( {k,\xi_{y},n_{\theta}} \right)}}}}}}}{\int{{\mathbb{d}^{3}r}\;{{\eta(r)}\left\lbrack {\frac{\alpha^{2}}{k_{0}} + {{\mathbb{i}}\left( {\rho - z_{0}} \right)}} \right\rbrack}^{- 1}{\int_{\pi}^{\pi}{{\mathbb{d}{{\phi exp}\left( {{\mathbb{i}}\;\phi\; n_{0}} \right)}}\exp\left\{ {{- {\mathbb{i}}}\;{r \cdot {{R(\phi)}\left\lbrack {{\hat{x}\sqrt{\left( {2k} \right)^{2} - \xi_{y}^{2}}} + {\hat{y}\xi_{y}}} \right\rbrack}}} \right\}}}}}} & (59)\end{matrix}$Where {tilde over (K)}(k,ξ_(p),n_(θ)) is given explicitly by

$\begin{matrix}{{{\overset{\sim}{K}\left( {k,\xi_{p},n_{\theta}} \right)} = {\frac{\mathbb{i}}{2}{A(k)}{\int_{0}^{\pi}{{\mathbb{d}\theta}\;{\exp\left( {{\mathbb{i}}\;\theta\; n_{\theta}} \right)}{\exp\left\lbrack {{- \frac{{\left( {2k} \right)^{2}\cos^{2}\theta} + {\xi_{p}^{2}\sin^{2}\theta}}{2}}\frac{\alpha^{2}}{2k^{2}}} \right\rbrack}}}}}{{\exp\left\lbrack {{- {\mathbb{i}}}\; z_{0}\sin\;\theta\sqrt{\left( {2k} \right)^{2} - \xi_{p}^{2}}} \right\rbrack}\sqrt{\left( {2k} \right)^{2} - \xi_{p}^{2}}\sin\;{\theta.}}} & (60)\end{matrix}$

Eq. (59) can be rewritten to use a Fredholm-type kernel κ(k, p, θ, r)such thatS(k,p,θ)=κη=∫d ³ rκ(k,p,θ,r)η(r).  (61)

Although Eq. (61) may not be strictly solvable, a least-squares solutionη⁺ can be found by minimizing a functional:

$\begin{matrix}{\eta^{+} = {{\underset{\eta}{\arg\;\min}{{S - {\kappa\;\eta}}}^{2}} = {\underset{\eta}{\arg\;\min}{\int_{0}^{\infty}{{\mathbb{d}k}{\int_{- \infty}^{\infty}{{\mathbb{d}p}{\int_{- \pi}^{\pi}{{\mathbb{d}\theta}{{{{S\left( {k,p,\theta} \right)} - {\kappa\;{\eta(r)}}}}^{2}.}}}}}}}}}} & (62)\end{matrix}$

The least-squares solution is then given by the pseudo-inverseη⁺=(κ^(†)κ)⁻¹κ^(†)S. While this solution is formally correct, theinverse (κ^(†)κ)⁻¹ can be difficult to compute in practice. Instead, wefind the least-squared error solution for the weighted Fourier transform{tilde over ({tilde over (η)}⁺ that, while not directly minimizing theerror relative to the measurements, still constrains the estimatedobject structure to be consistent with the data:

$\begin{matrix}{{\overset{\overset{\sim}{\sim}}{\eta}}^{+} = {{{\underset{\overset{\overset{\sim}{\sim}}{\eta}}{\arg\;\min}{{S - {K\;\overset{\overset{\sim}{\sim}}{\eta}}}}^{2}} + {\lambda{\overset{\overset{\sim}{\sim}}{\eta}}^{2}}} = {{\underset{\overset{\overset{\sim}{\sim}}{\eta}}{\arg\;\min}{\int_{0}^{\infty}{{\mathbb{d}k}{\int_{- \infty}^{\infty}{{\mathbb{d}\xi_{p}}{\sum\limits_{n_{\theta} = {- \infty}}^{\infty}{{{\overset{\overset{\sim}{\sim}}{S}\left( {k,\xi_{p},n_{\theta}} \right)} - {{K\left( {k,\xi_{p},n_{\theta}} \right)}{\overset{\overset{\sim}{\sim}}{\eta}\left( {k,\xi_{p},n_{\theta}} \right)}}}}^{2}}}}}}} + {\lambda{{{\overset{\overset{\sim}{\sim}}{\eta}\left( {k,\xi_{p},n_{\theta}} \right)}}^{2}.}}}}} & (63)\end{matrix}$

This least-squares solution keeps the object estimate consistent withEq. (57). Also included is a Tikhonov regularization term to stabilizethe solution, with regularization parameter λ. The solution {tilde over({tilde over (η)}⁺ is:

$\begin{matrix}{{{\overset{\overset{\sim}{\sim}}{\eta}}^{+}\left( {k,\xi_{p},n_{\theta}} \right)} = {\frac{{\overset{\overset{\sim}{\sim}}{S}\left( {k,\xi_{p},n_{\theta}} \right)}{K^{*}\left( {k,\xi_{p},n_{\theta}} \right)}}{{{K\left( {k,\xi_{p},n_{\theta}} \right)}}^{2} + \lambda}.}} & (64)\end{matrix}$

This least squares solution is a numerically simpler method ofestimating the object structure. Starting with data given by S(k, p, θ),we can compute {tilde over ({tilde over (S)}(k, ξ_(p), n_(θ)) using Eq.(52). Using Eq. (64) one can compute {tilde over ({tilde over(η)}⁺(k,ξ_(p), n_(θ)). Then Eq. (56) can be solved for {tilde over (η)}⁺by taking the discrete inverse Fourier transform of {tilde over ({tildeover (η)}⁺ with respect to n_(θ). Finally, a 3-D inverse Fouriertransform computes η⁺(r) from {tilde over (η)}⁺. In the limit that λ→0and all data are continuously available, this approach yields an exactsolution for η(r). In the more realistic scenario that a regularizedsolution is employed, a stable solution is obtained.

Simulation of Azimuthally-Scanned ISAM

As in the full-field OCT embodiment discussed above, theazimuthally-scanned algorithmic embodiment is now demonstrated byimplementing a simulation of the forward and inverse scattering in theradial OCT geometry. A synthetic object was created comprised ofpointlike scatterers. The simulated OCT data were calculated from theexact forward problem using Eq. (45), and then the regularized solutionof the inverse scattering solution was calculated using Eq. (64). Thesimulation is of a pseudo-three-dimensional object that is invariantalong the y-axis, so that the object is effectively two-dimensional.

The simulation was performed with lengths in units of the centralwavelength of the illumination. Typical center wavelengths for OCTimaging are between 800 nm and 1400 nm. The cross-section of thesimulated object was taken to be approximately 135 by 135 wavelengths.The Gaussian illumination beam was assumed to be focused 45 wavelengthsfrom the origin, with a beam waist width of 2.5 wavelengths. Thescatterers were placed 15 to 60 wavelengths from the origin at randomlyselected angles relative to the x-axis. The simulated source was takento have a Gaussian spectrum with a full-width-half-maximum (FVHM)fractional bandwidth of approximately 25%. Each of the scatterers hadthe same susceptibility.

The forward scattering problem was implemented by directly summing thecontribution of each point scatterer individually using Eq. (45). Thiswas achieved by summing for each sampled data point {tilde over(S)}(k,θ) the total collected backscattered amplitude for all of thescatterers at their respective positions r′ the amplitude f²(r′;k) asspecified in Eq. (44). Note that in Eq. (44) the exact phase rather thanthe Fresnel quadratic approximation was used to more accurately computethe data for a high numerical aperture beam. To show the equivalent OCTimage, the data was inverse Fourier transformed along the k axis,yielding S(r,θ). The resulting S(r,θ) is displayed in a polar plot inFIG. 9( a).

The dotted circle in the diagram indicates the radius at which theGaussian beam is focused. Note that the images of points located closerto the origin than the focus (inside the circle) curve towards theorigin, and the points located further from the origin than the focuscurve (outside the circle) away from the origin, as would be expected.

The inverse scattering problem was implemented using the approximatesolution embodied in Eq. (64). The data is given as {tilde over(S)}(k,θ). To utilize Eq. (64), several Fourier transform steps wereneeded. The inverse scattering algorithm was implemented using thefollowing steps:

-   -   The data {tilde over (S)}(k,θ) was Fourier transformed with        respect to 0 to yield {tilde over ({tilde over (S)}(k,n_(θ)).    -   The function {tilde over (K)}(k,θ) was calculated (using        ξ_(p)=0) and then Fourier transformed with respect to θ to yield        k(k, n_(θ)) as per Eq. (61).    -   The regularized {tilde over ({tilde over (η)}⁺(k,n_(θ)) was        calculated using Eq. 65.    -   {tilde over ({tilde over (η)}⁺(k,n_(θ)) was inverse Fourier        transformed with respect to n_(θ) to yield {tilde over        (η)}⁺(k,θ).    -   The {tilde over (η)}⁺(k,θ) was inverse Fourier transformed with        respect to k to yield η_(R) ⁺(l,θ), the Radon transform of        η⁺(x,z).    -   The inverse Radon transform of η_(R) ⁺(l,θ) was performed to        yield η⁺(x,z), the Tikhonov-regularized inverse.

The inverse Radon transform was used as a convenient way to convert fromthe polar representation of the Fourier transform {tilde over (η)}⁺(k,θ)to its inverse Fourier transform Cartesian counterpart η⁺(x,z), usingthe Fourier projection slice theorem. Unfortunately, manyimplementations of the inverse Radon transform, such as thefiltered-backprojection method that was used for this simulation, areslow, and therefore care will need to be exercised to ensure that thecomputational burden is not too great. Methods exist to implement theinverse Radon transform in O(N² log N) time, rather than the O(N³)typical of most filtered-backprojection inverse Radon transform methods.

The results of the inverse scattering computation are shown in FIG. 9(b). As can be seen, the blurred arcs corresponding to the point sourcesin the uncorrected OCT data are corrected to be pointlike when inversescattering is performed on the data. The algorithm correctly compensatesfor the range-dependent blur and curvature of the backscattered signal.Unlike in the translationally-scanned Gaussian beam or the full-fieldcases, the reconstructed image does not exhibit uniform resolution. Theresolution of the reconstruction depends on the distance from theorigin. Because the beam width is wide near the origin, points nearerthe origin than the focus are overlapped by the beam for many angles θ,so that the resolution of points near the origin is high. At the focus,the beam width is narrow and so points near the focus are also resolvedwell. Far from the origin, the beam is wide and points are onlyoverlapped by the beam for a narrow range of angles given by thedivergence of the beam, so that the resolution degrades with distancefrom the origin. Generally, the resolution is nearly constant betweenthe origin and the focus radius, and slowly degrades to a constantangular resolution at radii further than the focus. Therefore, the mostuseful resolution will be achieved for distances at or closer than thefocus radius.

To explore the range-dependent resolution further, a simulation of pointscatterers reconstructed with beams with various widths and focus radiiis now described with reference to FIG. 10. FIG. 10 has four parts, eachof which is the simulated resolution of point scatterers for beams ofdifferent widths. The marking next to each curve is the focus radius foreach simulated beam. The resolution is measured as the FWHM of thereconstructed point in the angular direction. Each graph relates theFWHM resolution to the distance from the axis to the simulated point.For small beam waists, as in parts (a), (b), and (c), the resolution isapproximately constant for radii closer than the focus radius. Furtherthan the focus the FVHM resolution starts to increase. For the widerbeams, the transverse resolution near the origin can be somewhat betterthan the width of the beam waist.

To examine the validity of the approximation made in Eq. 47 of smallfractional bandwidth, we simulate the reconstruction of point scatterersimaged with various source bandwidths. The simulated focus radius is 25wavelengths, and the beam widths are 1.5 and 5 wavelengths. FIG. 11shows the FVHM axial resolution as a function of fractional bandwidth.The resolution should be approximately half the reciprocal of thefractional bandwidth, to which the simulation conforms.

Phase Stability in ISAM

The increased resolution gained by the ISAM solution relies upon phasestable measurements. The phases in cross-correlation signals correspondto the measured relative delay between the reference and sample signalat particular temporal frequencies of the illumination. Theaforementioned methods rely on the phases of the cross-correlationmeasurements for a particular data set to correspond, i.e. the relativedelays measured for various positions of the illumination beam need toreflect the true delays for the beam to scatter off of layers inside theobject. Unfortunately, because of motions and changes in the objectduring the data acquisition, and thermal fluctuations that cause objectdimensions to change, these delays can vary during the data acquisitioninterval. This is a problem for almost all computed imaging modalities,but can be exacerbated for optical imaging because the small size of theillumination wavelength (often less than a micron) means that very smallmotions can cause relatively large fluctuations in phase. If the dataacquisition is brief, and the reference delay mechanism can produce areliably repeatable delay to an accuracy of less than a wavelength, thenphase stability between measurements can be easily achieved. Forexample, with the use of spectral detection for OCT, spectral-domain OCT(SD-OCT) seen in FIG. 6, we can be assured of phase stability withineach axial data set because the reference mirror is typically fixed andthe data acquisition is very rapid, typically occurring in fractions ofa second. Specifically, each axial acquisition is determined directlyfrom Fourier transform of the ensemble of spectral intensitymeasurements over the duration of the exposure time on a CCD camera.Thus, relative phases between adjacent reflections in the sample arefixed relative to each other and the reference for a single axialacquisition. Furthermore, if adjacent axial scans may be captured fastenough to avoid some minimum amount of phase drift then an accuratereconstruction is possible. Phase drift can occur in a system formultiple reasons including thermal changes, galvanometer or stagepositioning accuracy, and system or sample jitter. The greater the timeinterval between scans, the more likely it is that random phase errorswill occur. Adjacent axial scans in a single cross-sectional scan arethus less likely to suffer from distortions due to random phase jitterthan adjacent axial scans from multiple cross-sectional scans.

Object reconstruction requires the phase to be stable relative to allaxial scans of a 3D acquisition. There are several ways to achieve phasestability, whether it is one of many hardware and environmentalconfigurations, implementations of reference objects, or algorithmicdevelopments based on priors. In conventional scanned beam OCT, it hasbeen shown by Ralston et al., Phase Stability Technique for InverseScattering in Optical Coherence Tomography, IEEE International Symposiumon Biomedical Imaging., pp. 578-581, (2006), incorporated herein byreference, that one such possible method to achieve 3D phase stabilityin OCT for reconstruction of the inverse scattering solution is to use aflat reference reflector such as a microscope coverslip. Because thecoverslip typically remains in contact with the object, its positionrelative to the object is fixed, and therefore can serve to indicate adelay that is consistent between axial scans. Such a method offersadvantages over expensive environmental controllers and extremely fastacquisition hardware. Further, we develop an algorithm for tracking theair-glass interface. Other interfaces that are fixed relative to theobject being measured can also be used, such as the interior or exteriorsurface of the transparent wall of a catheter.

The acquired SD-OCT signal can be represented after dispersioncorrection as a function of transverse position and wave number,S(r₀,k), where the wave numbers k are related to the frequencies ω bythe dispersion relation k(ω)=ω n/c, and n is the index of refraction.

We present a method that analyzes each axial scan individually andapplies a phase to compensate variations of the position of the samplerelative to the illumination beam. We place a reflective surface such asa microscope coverslip on top of the sample to act as a referencesurface, which is used to infer the delay to the top surface of thesample. The signal for an arbitrary axial line of data can berepresented as S(k), a function of k. We assume that there is a range ofdistances along the illumination beam z_(min) to z_(max) for which thesignal reflected from the coverslip is known to reside in every axialscan. The inverse Fourier transform of S(k) is computed as S_(c)(z), andthe signal corresponding to the reflection is contained in the samplesS_(c)(z) for z_(min)≦z≦z_(max). The spatial spectrum of the reflectionis computed as the Fourier transform of S_(c)(z) over the windowz_(min)≦z≦z_(max), which is called {tilde over (S)}_(c)(k).

Because the signal contained in {tilde over (S)}_(c)(k) corresponds to asingle reflection, it is reasonable to model it as {tilde over(S)}_(c)(k)=A(k)e^(eφ(k)), where the phase function φ(k)=φ₀+kd, where φ₀is an arbitrary phase and d is the true position of the surface wherethe reference reflection occurs. Because of the motion of the sample,the actual phase arg {tilde over (S)}_(c)(k)=φ′(k). By multiplying theaxial scan data S(k) by the correction factor e^(i[φ(k)−φ′(k)]), thephase of the axial scan can be adjusted to place the reflection at itstrue known position d.

We model the phase φ′(k), as a Taylor series around a center frequencyk₀:

${{\phi^{\prime}(k)} = {{{\phi^{\prime}\left( k_{0} \right)} + {\left( {k - k_{0}} \right)\frac{\partial\phi^{\prime}}{\partial k}}}❘_{k = k_{0}}{+ \ldots}}}\mspace{11mu},$To utilize this model, we must estimate the value of ∂φ′/∂k|_(k=k) ₀from the data function φ′(k). The function φ′(k) is wrapped to the range−π to π, so calculating the derivative directly from the wrapped datawill incorrectly incorporate the 2π jumps into the estimate. Instead, weform the unwrapped φ_(w)(k) by removing 2π discontinuities from φ′(k).The estimate then becomes

$\frac{\partial\phi^{\prime}}{\partial k}❘_{k = k_{0}}{\approx \frac{{\phi_{w}\left( k_{2} \right)} - {\phi_{w}\left( k_{1} \right)}}{k_{2} - k_{1}}}$where k₁<k₀<k₂, with the frequencies k₁ and k₂ chosen to span theillumination spectrum, typically with k₁ and k₂ corresponding to thefrequencies at which the power spectral density is half of that at thepeak.

Once φ′(k₀) and ∂φ′/∂k|_(k=k) ₀ are known, the empirical φ′(k) can becomputed, and the corrected axial scan spectrumS′(k)=S(k)e^(i[φ(k)−φ(k′)]). This corrected axial scan data will bemodified such that the position of the reference reflection is always atthe same location on the axial scan, removing the effective longitudinalrelative motion between the sample and the scanned beam. For this methodto work properly, the reference object must be located for each axialscan, otherwise that axial scan could contribute to a poorreconstruction. Furthermore, refinements to this method could utilizehigher order terms of the series for φ′(k), which would account forinstrument dispersion as well as motion.

Experimental Examples of ISAM Imaging

A tissue phantom consisting of a collection of titanium dioxidescatterers having a mean diameter of 1 μm and uniformly suspended insilicone was imaged using low-coherence interferometry and a 0.05 NAobjective. FIG. 12 displays cross-sections through an unprocessed dataset (left) and ISAM reconstruction (right) of a volume 360 μm×360 μm(transverse)×2000 μm (axial). FIG. 12 a-f contain three pairs of en facesections for both the unprocessed data (FIG. 12 a-c) and the ISAMreconstructions (FIG. 12 d-j). The distances from the en face sectionplanes to the focus, located at z=0, are z=1100 μm (FIGS. 12 a and 12a), z=475 μm (FIGS. 12 b and 12 e), and z=−240 μm (FIGS. 12 c and 12 f).These sections show that the reconstruction has resolved the scatterersthroughout a range of depths over nine times the Rayleigh range from thefocus, where the Rayleigh range is commonly defined as half of thedepth-of-field, or what is considered in-focus in optical imagingsystems. In the unprocessed data, the interference between the signalsscattered from adjacent scatterers is evident. Our method properlyaccounts for the diffraction of the beam, and so separates thesuperimposed signals from the scatterers to form well-resolved pointimages on all planes.

Human tumor tissue was resected and imaged ex vivo. Sections were markedwith India ink after imaging and before embedding to register locations.FIG. 13 includes en face planes (Slices A and B) of the unprocessed data(FIG. 13 a) where the beam diffraction effects are evident, the computedISAM reconstructions (FIG. 13 b), and images of corresponding registeredhistological sections (FIG. 13 c). Although embedding, sectioning, andstaining of tissue can disrupt features to some degree, the registeredhistological sections provide prominent features for comparison. In thereconstruction, boundaries between adipose cells and the margin betweenadipose and fibrous tissue are clearly identified, with a strongcorrespondence to histology. While the histological images were obtainedby destroying the sample, ISAM could readily be applied for in vivoapplications because signal collection is in the back-scatteredepi-direction.

Real-Time Cross-Sectional Processing

In order to provide the benefits of spatially-invariant resolutionattainable by means of the hitherto described features of ISAM and toobtain immediate feedback in time-critical situations or for monitoringtransient dynamics, methods are now described for real-time realizationof ISAM computations.

As described in detail above, an object being imaged by spectral-domainOCT (SD-OCT) has a susceptibility represented by η(r_(∥),z), where r_(∥)is the transverse position and z is the longitudinal position. Thecollected SD-OCT signal is represented by {tilde over (S)}(r_(∥),ω),whose arguments are the transverse position of the beam r_(∥) and thefrequency ω. After correcting for dispersion, the signal is representedby {tilde over (S)}(r_(∥),k), where k is the uniformly spacedwavenumber. The Fourier transform of {tilde over (S)}(r_(∥),k) withrespect to k is S(r_(∥),t). Introducing a coordinate change from t to t′such that t′=0 coincides with the delay of the focal plane results in asignal S(r_(∥), t′). The 3-D Fourier transform of S(r_(∥),t′) is {tildeover ({tilde over (S)}(Q,k).

The algorithm for ISAM is constructed as described above. A solution forhigh-numerical aperture optics is implemented, namely, a relation isderived that associates the Fourier transform of the object η with theFourier transform of the signal S. The 2-D Fourier transform of thespectral-domain OCT signal {tilde over (S)}(r_(∥),k), is given by theexpression

$\begin{matrix}{{{\overset{\overset{\sim}{\sim}}{S}\left( {Q,k} \right)} = {\frac{k^{2}}{\alpha^{2}}i\; 2\pi^{2}{A(k)}\frac{e^{{- 2}{{ik}_{z}{({Q/2})}}z_{0}}}{k_{z}\left( {Q/2} \right)}e^{\frac{{- u^{2}}Q^{2}}{2k^{2}}}{\overset{\overset{\sim}{\sim}}{\eta}\left\lbrack {Q,{2{k_{z}\left( {Q/2} \right)}}} \right\rbrack}}},} & (65)\end{matrix}$where {tilde over ({tilde over (η)} is the 3-D Fourier transform of η,the argument Q is the transverse wavevector, k is the wavenumber,k_(z)(q)=√{square root over (k²−q²)}, z₀ is the fixed position of thebeam focus, A(k) is the square root of the power spectral density,α=π/NA, and NA is the numerical aperture of the beam. The correspondingTikhonov regularized solution, {tilde over ({tilde over (η)}⁺, takes theform

$\begin{matrix}{{{{\overset{\sim}{\overset{\sim}{\eta}}}^{+}\left( {Q,\beta} \right)} = \left\lbrack \frac{{f^{*}\left( {Q,k,\beta} \right)}{\overset{\sim}{\overset{\sim}{S}}\left( {Q,k} \right)}}{{{f\left( {Q,k,\beta} \right)}}^{2} + {2\lambda\;{k/{k_{z}\left( {Q/2} \right)}}}} \right\rbrack_{k - {\frac{1}{2}\sqrt{\beta^{2} - Q^{2}}}}}{where}} & (66) \\{{{f\left( {Q,k,\beta} \right)} - {\frac{k^{2}}{\alpha^{2}}i\; 2\;\pi^{2}{A(k)}\frac{e^{{- 2}{{ik}_{z}{({Q/2})}}z_{0}}}{k_{z}\left( {Q/2} \right)}e^{\frac{\alpha^{2}Q^{2}}{4k^{2}}}}},} & (67)\end{matrix}$β is the longitudinal frequency coordinates of the object, and λ is theregularization constant. Rearranging the terms of the Fourier space intomultiplicative factors of magnitude, {tilde over ({tilde over (B)}(Q,k),and phase, e^(i(2k) ^(z) ^((Q/2)z) ⁰ ^(+x/2)), the pseudo inverse can berewritten as

$\begin{matrix}{{{{\overset{\overset{\sim}{\sim}}{\eta}}^{+}\left( {Q,\beta} \right)} = {{{\overset{\overset{\sim}{\sim}}{B}\left( {Q,k} \right)}e^{i{({{2{k_{z}{({Q/2})}}z_{0}} + {\pi/2}})}}{\overset{\overset{\sim}{\sim}}{S}\left( {Q,k} \right)}}❘_{k = {\frac{1}{2}\sqrt{\beta^{2} + Q^{2}}}}}},{where}} & (68) \\{{\overset{\overset{\sim}{\sim}}{B}\left( {Q,k} \right)} = {\frac{{- \frac{k}{\alpha^{2}}}\pi^{2}{A(k)}e^{\frac{\alpha^{2}Q^{2}}{4k^{2}}}}{{\frac{k^{3}}{\alpha^{4}}\frac{2\pi^{4}{A^{2}(k)}}{k_{z}\left( {Q/2} \right)}e^{\frac{{- \alpha^{2}}Q^{2}}{2k^{2}}}} + \lambda}.}} & (69)\end{matrix}$Without loss of generality, the same origin is designated for the depthcoordinates z as the time delay coordinates t. The t=0 delay correspondsto the z=0 plane and coincide at the path length matched delay of thereference. Additionally, for a coordinate change, the t′=0 delaycorresponds to the z′=0 plane and coincide at the focal plane, wherethen z₀ will equal zero, and equation (66) reduces to

$\begin{matrix}{{{{\overset{\sim}{\overset{\sim}{\eta}}}^{+}}^{\prime}\left( {Q,\beta} \right)} = {{{\overset{\overset{\sim}{\sim}}{B}\left( {Q,k} \right)}e^{i\;{\pi/2}}{\overset{\overset{\sim}{\sim}}{S}\left( {Q,k} \right)}}❘_{k = {\frac{1}{2}\sqrt{\beta^{2} + Q^{2}}}}.}} & (70)\end{matrix}$

A coordinate change from t to t′ is achieved by circularly shifting thedata. The origin of the t coordinates is at the time corresponding tothe arrival of the reference light. The origin of the t′ coordinates isshifted such that t′=0 is at the focal plane. A zero pad of 2|t₀| rowsprevents aliasing of the shifted data, where t₀ is the number of delaysamples from the focus to t′=0. The model in FIG. 15( a) and the grid inFIGS. 15( b) and (c) describe visually in 2D the ISAM resampling seen inequation (70).

A shift of t₀ in S(r_(∥),t) is used to make the t=0 delay coincide withthe z=0 plane as will be seen in the algorithm. S(r_(∥),t+t₀) has theFourier transform relation {tilde over ({tilde over (S)}(Q,k)e^(−kQt) ⁰. If we substitute this into the equation.

$\begin{matrix}{{\overset{\overset{\sim}{\sim}}{\eta}\left( {Q,\beta} \right)} = {{{\overset{\sim}{\overset{\sim}{B}}\left( {Q,k} \right)}e^{2{{ik}_{z}{({Q/2})}}z_{0}}{\overset{\overset{\sim}{\sim}}{S}\left( {Q,k} \right)}e^{- {kQt}_{0}}}❘_{k - {\frac{1}{2}\sqrt{\beta^{2} + Q^{2}}}}}} & (71) \\{{{\overset{\overset{\sim}{\sim}}{\eta}}^{+}\left( {Q,\beta} \right)} = {{{\overset{\overset{\sim}{\sim}}{B}\left( {Q,k} \right)}{ie}^{i({{2{k_{z}{({Q/2})}}z_{0}} - {kt}_{0}}\}}{\overset{\overset{\sim}{\sim}}{S}\left( {Q,k} \right)}}❘_{k = {\frac{1}{2}\sqrt{\beta^{2} + Q^{2}}}}}} & (72)\end{matrix}$This equation may be simplified by having the relative time offset t₀and the focus offset z₀ coincide at t₀=0 and z₀=0.

FIG. 16 shows a data flow diagram for the prototype ISAM algorithm for2D imaging where in the paraxial limit r_(∥) now describes a singletransverse dimension. It is noted that the other transverse dimensionmay be solved separately in so much as the paraxial approximation holds.The prototype algorithm was designed in Matlab, a mathematical softwarepackage, where all the calculations were done in double precision,64-bit. The digitized spectra S_(rω) is a function of the M discretepositions r in the plane perpendicular to the beam axis and the Ndiscrete frequencies ω. A non-uniform interpolation is needed to map thesampled frequencies ω to the sampled wavenumbers k. Similarly, the ISAMsolution requires a non-uniform interpolation mapping wavenumbers k tolongitudinal frequency coordinates of the object β. The cubic B-splineis chosen as the non-uniform resampling interpolator; although, awindowed sinc interpolator, such as the prolate-spheroidal interpolator,may be chosen to select specific convergence characteristics.Unfortunately, the frequency response for many non-uniform interpolationprocedures drops in performance for frequencies higher than half theNyquist rate. To mitigate this effect, each spectrum is interpolated byperforming a fast Fourier transform (FFT), padding with N zeros, andperforming an inverse FFT (IFFT). The interpolated spectra are stored ina buffer with 2N rows and M columns. Thus, the new interpolated,digitized spectra S_(rω′) is a function of the sampled positions r andthe new sampled frequencies ω′. Similarly, S_(rk) is interpolated by afactor of 2 to get S_(rk′) as a function of r and the new uniformlysampled wavenumbers k′.

The non-uniformly sampled spectrum in spectral-domain OCT can becorrected in a fashion similar to material dispersion correction byresampling the Fourier spectrum from ω′ to k. A desired reindexing arrayi_(n) is based primarily on calibrated, second- and third-ordercorrection factors α₂ and α₃, respectively.

$\begin{matrix}{{i_{n} = {{2n} + {\alpha_{2}\left( {\frac{2n}{N} - \omega_{ctr}} \right)}^{2} + {\alpha_{3}\left( {\frac{2n}{N} - \omega_{ctr}} \right)}^{3}}},} & (73)\end{matrix}$where n is an integer between 0 and N−1, 2N is the number of samples ofthe interpolated spectrum S_(rω′), and ω_(ctr) is the calculatedcentroid of the spectrum on a scale from 0 to 1. α₂ and α₃ are setthrough a calibration routine. A mirror models a perfect reflector witha minimized width of the point spread function (PSF). Thus, α₂ and α₃are adjusted such that the imaged PSF is indeed minimized.

FIG. 16 shows a data flow chart of the prototype algorithm which hasbeen broken down into the eight steps listed below.

Step 1 The spline coefficients and interpolations are computed asdescribed above. The result is stored in an array S_(rk′) where k′ isthe uniformly sampled wavenumber.

Step 2 The 1-D FFT of the columns of S_(rk′) is taken to get the signalS_(rt) where t is the sampled time delay of the optical signal. TheHermitian symmetric values, where t<0, are removed. There will be noambiguity, if the reference arm is placed such that the reference lightpulse arrives at the spectrometer before the sample light pulses. Thecomplex analytic OCT signal is represented by S_(rt).Step 3 A phase correction procedure is implemented to ensure phasestability for the ISAM reconstruction. A glass coverslip can be placedon top of the sample to track and reduce system instabilities. The phasecorrected signal is represented as Ŝ_(rt).Step 4 The contribution of the constant spectral intensity is removed bysetting Ŝ_(rt) to zero, for values near t=0. Then, a coordinate changefrom t to t′ is made by circularly shifting the rows such that the focalplane lies at t′=0, which coincides with z′=0 where the solution isderived in equation (72). The data Ŝ_(rt′) is padded with 2|t₀| rows ofzeros to prevent aliasing of the shifted data, where t₀ is the number oftime delay samples from the focus to the center of Ŝ_(rt).Step 5 The 2-D IFFT of Ŝ_(rt′) is taken to get S_(qk).Step 6 The 2-D ISAM resampling grid, spline coefficients, and theinterpolations are calculated. Then the result is multiplied by {tildeover ({tilde over (B)}(Q,k)e^(iπ/2) to get η_(qβ).Step 7 The 2-D FFT of η_(qβ) is taken to get η_(rz′), where z′=0 at thefocal plane.Step 8 A coordinate change from z′ to z is made by circularly shiftingthe rows such that the reference delay lies at the top of the image andcorresponds to the z=0 plane. Then, the magnitude of the ISAMreconstruction η_(rz) is calculated resulting in ISAM amplitude image|η_(rz)|.

There are a number of operations hindering the performance of theprototype ISAM algorithm. Using C++ instead of Matlab allows for anumber of speed improvements. The 64-bit operations such as the FFT andinterpolations can be replaced with 32-bit operations with a small, buttolerable, increase of quantization noise. A speed advantage is gainedby replacing the ‘for’ loops within Matlab scripts by vectorized IntelSSE (Streaming SIMD Extentions) commands and/or compiled ‘for’ loops.Time-expensive, built-in Matlab interpolation and FFT functions arereplaced with hardware optimized functions as found in the Intel MathKernel Library (MKL). An FFT of the real spectrum is implemented usingthe customized real-to-complex FFT in the MKL. The resampling stepcorrects the depth dependent defocus and is crucial for the performanceof the algorithm. Although, the filter in equation (72), {tilde over({tilde over (B)}(Q,k)e^(iπ/2), can be highly singular which introducesnoise, hence the need for regularization. Furthermore, applying thefilter provides a quantitative, but not a significant qualitative changeto the data. Thus, equation (72) is reduced to an unfiltered solution

$\begin{matrix}{{{\overset{\overset{\sim}{\sim}}{\eta}}^{+^{\prime\prime}}\left( {Q,\beta} \right)} = {{\overset{\overset{\sim}{\sim}}{S}\left( {Q,k} \right)}❘_{k = {\frac{1}{2}\sqrt{\beta^{2} + Q^{2}}}}}} & (74)\end{matrix}$without degradation of the Fourier space resampling. The focal plane isfixed such that t₀=0 by imposing a restriction that the focus be placedat the center of the image. Therefore, the complex analytic signal doesnot need to be padded with any zeros, and thus the computation time ofthe FFT is reduced because the FFT is always a power of two. While theprototype functions utilized dynamically allocated memory, the real-timeprogram allocates memory prior to imaging time. Similarly, a table ofresampling coefficients for the dispersion compensation and the ISAMreconstruction are pre-calculated and stored in memory. The prototypealgorithm interpolated the data by two to mitigate the high frequencyroll-off of the non-uniform interpolator. Although the interpolation hasgood frequency response, it is not necessary, especially when speed isan issue. The phase stabilization procedures, which might be needed forlong acquisition periods, are not necessary for 2-D imaging since thisSD-OCT system provides good phase stability over the short duration ofthe scan, about 17 ms.

The key computation for ISAM is the resampling in the Fourier space. Thenew design is an efficient routine which requires two interpolations ofthe columns, one one-dimensional (1D) real-to-complex (R2C) fast Fouriertransform (FFT) of the columns, and two 2D FFTs. The FFT routine fromthe Intel Math Kernel Library (MKL) was used for all 1D and 2Dtransforms. The ID 2048-point R2C FFT is calculated for every column ofthe 512×2048 real-valued data, while the 2D FFT and 2D IFFT arecalculated for the 512×1024 complex values.

The reindexing array i_(n) and the corresponding cubic B-splineinterpolation table is pre-computed and stored in random access memory(RAM) for repeated use. The coefficients for the cubic B-spline arepre-computed. The integer part of the index used for calculation of thein the cubic B-spline is given by

$\begin{matrix}{{a_{x}\lbrack n\rbrack} = \left\{ {\begin{matrix}{{\left\lfloor i_{n} \right\rfloor + x},} & {0 \leq {\left\lfloor i_{n} \right\rfloor + x} \leq {N - 1}} \\{0,} & {{\left\lfloor i_{n} \right\rfloor + x} < 0} \\{{N - 1},} & {{\left\lfloor i_{n} \right\rfloor + x} > {N - 1}}\end{matrix},{x = {- 1}},0,1,{{2\mspace{14mu}{and}\mspace{14mu} 0} \leq n < N}} \right.} & (75)\end{matrix}$The fractional coefficients are given byf _(n) =i _(n) −└i _(n)┘, 0≦n<N  (76)andb ⁻¹ [n]=(1−f _(n))³/6b ₀ [n]=(4−6f _(n) ²+3f _(n) ³)/6,b ₁ [n]=(1+3f _(n)+3f _(n) ²−3f _(n) ³)/6, 0≦n<N  (77)b ₂ [n]=f _(n) ³/6

Next, the Fourier resampling indices of ISAM are pre-calculated. Theconstants which specify the axial and transverse bandwidths of theobject, based on system parameters, are β_(min), β_(max), q_(min), andq_(max), respectively. These constants are selected in accordance to thespecific bandwidth parameters of the system and describe the boundariesof the Fourier space shown in FIG. 15( c). By defining the longitudinalbandwidth parameters of the object rather than the wavenumber, we canavoid computationally costly zeropadding of the resampled solution.However, a small loss of spectral information across the β_(min)boundary may reduce the signal-to-noise ratio, but only marginally. Inthis case, β_(min) and β_(max) are set to be the boundaries of theoptical bandwidth, q_(min) is set to zero, and q_(max) is set to themaximum transverse scanning frequency. More important than the exactvalues for β_(min), β_(max), q_(min), and q_(max), is the ratio of thecorresponding transverse and longitudinal bandwidths. The maximum andminimum wavenumber are calculated for the region of interest,k _(min)=β_(min)/2,  (78)k _(max)=0.5√{square root over (β_(max) ² +q _(max) ²)}.  (79)A range of values for q[m] and β[n] is created in the 2-D Fourier spaceand the resampling grid kq[m,n] is pre-calculated. Notice here that Mand N′=N/2 are assigned according to number of rows and columns in thecomplex analytic sampled Fourier data.

$\begin{matrix}{{q\lbrack m\rbrack} = \left\{ {\begin{matrix}{{m\frac{2q_{\max}}{M}},} & {0 < m \leq {M/2}} \\{{\left( {m - M} \right)\frac{2q_{\max}}{M}},} & {{M/2} < m \leq M}\end{matrix},} \right.} & (80) \\{{{\beta\lbrack n\rbrack} = {{{n\left( {\beta_{\max} - \beta_{\min}} \right)}/N^{\prime}} + \beta_{\min}}},{0 \leq n < N^{\prime}},} & (81) \\{{{{kq}\left\lbrack {m,n} \right\rbrack} = {{\frac{N^{\prime}}{k_{\max} - k_{\min}}\left( {{0.5\sqrt{{\beta\lbrack n\rbrack}^{2} + {q\lbrack m\rbrack}^{2}}} - k_{\min}} \right)} + 1}},{0 < n < N^{\prime}},{0 < m < M}} & (82)\end{matrix}$The kq[m,n] grid is used to calculate a table of values for the cubicB-spline coefficients. The values are calculated as shown above excepteach column of resampling parameters is different and therefore mustalso be saved. FIG. 1( b) shows the plot of the kq[m,n] grid where thecurved lines represent the constant values of β[n]. This 2D grid is usedto calculate another table of cubic B-spline coefficients. The splinevalues are calculated as shown above except each column of theresampling parameters is different and therefore 2D spline coefficientsare stored in memory.

Similar to the spline coefficient calculations shown above, thereindexing array kq[m,n] and the corresponding cubic B-splineinterpolation coefficient table is pre-computed and stored in randomaccess memory (RAM) for repeated use. The integer part of the index usedfor calculation of the in the cubic B-spline is given by

$\begin{matrix}{{a_{q,x}^{\prime}\left\lbrack {m,n} \right\rbrack} = \left\{ {{\begin{matrix}{{\left\lfloor {{kq}\left\lbrack {m,n} \right\rbrack} \right\rfloor + x},} & {0 \leq \left\lfloor {{{{kq}\left\lbrack {m,n} \right\rbrack} + x} \leq {N^{\prime} - 1}} \right\rfloor} \\{0,} & {{\left\lfloor {{kq}\left\lbrack {m,n} \right\rbrack} \right\rfloor + x} < 0} \\{{N^{\prime} - 1},} & {{\left\lfloor {{kq}\left\lbrack {m,n} \right\rbrack} \right\rfloor + x} > {N^{\prime} - 1}}\end{matrix};{x = {- 1}}},0,1,{2;{0 \leq n < N^{\prime}};{0 \leq m < M}}} \right.} & (83)\end{matrix}$The fractional coefficients are given byf _(m,n) =kq[m,n]−└kq[m,n]┘, 0≦n<N′; 0≦m<M  (84)andb′ _(q,−1) [m,n]=(1−f _(m,n))³/6b′ _(q,0) [m,n]=(4−6f _(m,n) ²+3f _(m,n) ³)/6b′ _(q,1) [m,n]=(1+3f _(m,n)+3f _(m,n) ²−3f _(m,n) ³)/6, 0≦n<N′;0≦m<M  (85)b′ _(q,2) [m,n]=f _(m,n) ³/6Using the pre-calculated tables, a flow diagram of the real-timealgorithm is shown in FIG. 17.

Here S_(rω)[m,n] is the raw interferometric data captured from thecamera and has M columns and N rows. In this implementation, M=512columns and N=2048 rows.

Step 1 The pre-calculated table is used to perform the interpolation asfollows.S _(rk) [m,n]=S _(rω) [m,a ⁻¹ {n}]b ⁻¹ {n}+S _(rω) [m,a ₀ {n}]b ₀ {n}+S_(rω) [m,a ₁ {n}]b ₁ {n}+S _(rω) [m,a ₂ {n}]b ₂ {n}  (86)for all integers 0≦n<N and 0≦m<M.Step 2 The real-to-complex 1-D FFT routine from the Intel Math KernelLibrary (MKL) is used on all the columns.

$\begin{matrix}{{{S_{rt}\left\lbrack {m,n} \right\rbrack} = {\sum\limits_{k = 0}^{N - 1}{{S_{rk}\left\lbrack {m,k} \right\rbrack}e^{{- \frac{2\pi\; i}{N}}{kn}}}}},{0 \leq {{nN}\mspace{14mu}{and}\mspace{14mu} 0} \leq m < M}} & (87)\end{matrix}$The real-to-complex FFT will compute N/2 complex values. The new numberof rows of the complex data is given by N′=N/2.Step 3 The contribution of the noise from the average spectral intensityon the detector is removed by setting S_(rt)[m,n] equal to zero at thet=0 plane. Also, S_(rt)[m,n] is circularly shifted by half such that thefocus will be the new t′=0 plane.

$\begin{matrix}{{S_{{rt}^{\prime}}\left\lbrack {m,n} \right\rbrack} = \left\{ {{\begin{matrix}{S_{rt}\left\lbrack {m,{n + {N^{\prime}/2}}} \right\rbrack} & {0 \leq n < {N^{\prime}/2}} \\0 & {{N^{\prime}/2} \leq n < {{N^{\prime}/2} + 2}} \\{S_{rt}\left\lbrack {m,{n - {N^{\prime}/2}}} \right\rbrack} & {{{N^{\prime}/2} + 2} \leq n < N^{\prime}}\end{matrix}\mspace{14mu}{and}\mspace{14mu} 0} \leq m < {M.}} \right.} & (88)\end{matrix}$Step 4 The complex 2-D inverse FFT (IFFT) of the complex analytic signalS_(rt′)[m,n] is calculated

$\begin{matrix}{{{S_{qk}\left\lbrack {m,n} \right\rbrack} = {\frac{1}{{MN}^{\prime}}{\sum\limits_{r = 0}^{M - 1}{\sum\limits_{t = 0}^{N^{\prime} - 1}{{S_{{rt}^{\prime}}\left\lbrack {r,t} \right\rbrack}e^{\frac{2\pi\; i}{N^{\prime}}{nt}}e^{\frac{2\pi\; i}{M}{mr}}}}}}},{0 \leq n < {N^{\prime}\mspace{14mu}{and}\mspace{14mu} 0} \leq m < {M.}}} & (89)\end{matrix}$Step 5 The pre-calculated table is used to perform the cubic B-splineinterpolation as follows.η_(qβ) [m,n]=S _(qk) [m,a′ _(q,−1) {m,n}]b _(q,−1) {m,n}+S _(qk) [m,a′_(q,0) {m,n}]b _(q,0) {m,n}+S _(qk) [m,a′ _(q,1) {m,n}]b′ _(q,1) {m,n}+S_(qk) [m,a′ _(q,2) {m,n}]b′ _(q,2) {m,n},0≦n<N′ and 0≦m<M,  (90)where the calculated cubic B-spline coefficients are from the lookuptable.Step 6 The complex 2-D FFT of the Fourier transformed object η_(qβ)[m,n]is calculated

$\begin{matrix}{{{\eta_{{rz}^{\prime}}\left\lbrack {m,n} \right\rbrack} = {\sum\limits_{q = 0}^{M - 1}{\sum\limits_{\beta = 0}^{N^{\prime} - 1}{{\eta_{q\;\beta}\left\lbrack {q,\beta} \right\rbrack}e^{{- \frac{2\pi\; i}{N^{\prime}}}\beta\; n}e^{{- \frac{2\pi\; i}{M}}{qm}}}}}},{0 \leq n < {N^{\prime}\mspace{14mu}{and}\mspace{14mu} 0} \leq m < {M.}}} & (91)\end{matrix}$Step η_(rz′)[m,n] is circularly shifted such that the focus is in themiddle of the image,

$\begin{matrix}{{\eta_{rz}\left\lbrack {m,n} \right\rbrack} = \left\{ {\begin{matrix}{\eta_{{rz}^{\prime}}\left\lbrack {m,{n + {N^{\prime}/2}}} \right\rbrack} & {0 \leq n < {N^{\prime}/2}} \\{\eta_{{rz}^{\prime}}\left\lbrack {m,{n - {N^{\prime}/2}}} \right\rbrack} & {{N^{\prime}/2} \leq n < N^{\prime}}\end{matrix},} \right.} & (92)\end{matrix}$then the magnitude |η_(rz)[m,n]| is displayed.

In various embodiments of the present invention, the disclosed methodsdetermining the three-dimensional susceptibility of a sample may beimplemented as a computer program product for use with a computersystem. Such implementations may include a series of computerinstructions fixed either on a tangible medium, such as a computerreadable medium (e.g., a diskette, CD-ROM, ROM, or fixed disk) ortransmittable to a computer system, via a modem or other interfacedevice, such as a communications adapter connected to a network over amedium. The medium may be either a tangible medium (e.g., optical oranalog communications lines) or a medium implemented with wirelesstechniques (e.g., microwave, infrared or other transmission techniques).The series of computer instructions embodies all or part of thefunctionality previously described herein with respect to the system.Those skilled in the art should appreciate that such computerinstructions can be written in a number of programming languages for usewith many computer architectures or operating systems. Furthermore, suchinstructions may be stored in any memory device, such as semiconductor,magnetic, optical or other memory devices, and may be transmitted usingany communications technology, such as optical, infrared, microwave, orother transmission technologies. It is expected that such a computerprogram product may be distributed as a removable medium withaccompanying printed or electronic documentation (e.g., shrink wrappedsoftware), preloaded with a computer system (e.g., on system ROM orfixed disk), or distributed from a server or electronic bulletin boardover the network (e.g. the Internet or World Wide Web). Of course, someembodiments of the invention may be implemented as a combination of bothsoftware (e.g., a computer program product) and hardware. Still otherembodiments of the invention are implemented as entirely hardware, orentirely software (e.g., a computer program product).

The embodiments of the invention heretofore described are intended to bemerely exemplary and numerous variations and modifications will beapparent to those skilled in the art. All such variations andmodifications are intended to be within the scope of the presentinvention as defined in any appended claims.

1. A method for determining the three-dimensional susceptibility of asample, the method comprising: a. providing a source of a beam ofradiation characterized by a temporally dependent spectrum and a localpropagation axis defined at a locus of incidence with respect to asurface of the sample; b. irradiating the surface in a planecharacterized by a fixed displacement along the propagation axis of thebeam; c. superposing scattered radiation from the sample with areference beam derived from the source of the beam of radiation toprovide an interference signal; d. deriving a forward scattering modelrelating measured data to structure of an object; and e. solving aninverse scattering problem based upon the forward scattering model andthe interference signal to infer a three-dimensional structure of thesample.
 2. A method in accordance with claim 1, wherein the step ofproviding a source of a beam of radiation includes providing a partiallycoherent source.
 3. A method in accordance with claim 1, wherein thestep of deriving a forward scattering model further includes relating atransform of measured data to a transform of structure of the sample. 4.A method in accordance with claim 3, further comprising a step oftransforming the interference sign into a transform of the interferencesignaling a transform domain.
 5. A method in accordance with claim 1,wherein the step of solving an inverse scattering problem is performedin substantially real time.
 6. A method in accordance with claim 1,wherein the step of irradiating includes focusing radiationsubstantially at the surface of the sample.
 7. A method in accordancewith claim 6, wherein the step of focusing includes illuminating thesample through a microscope objective.
 8. A method in accordance withclaim 1, wherein the step of irradiating the sample includes deliveringthe beam to the sample by at least one of a catheter, a needle, a probe,and endoscope, and a laparoscope.
 9. A method in accordance with claim1, further comprising sweeping the spectrum of the source of radiationas a function of time.
 10. A method in accordance with claim 1, whereinthe step of providing a source includes providing a tunable laser.
 11. Amethod in accordance with claim 1, wherein the step of superposingincludes configuring arms of the beam of radiation in an interferometer.12. A method in accordance with claim 11, wherein the step ofsuperposing includes configuring arms of the beam in one of a Michelson,a Mach-Zehnder, and a Fizeau interferometer.
 13. A method in accordancewith claim 1, wherein the step of superposing the scattered radiationand reference beam includes superposition with a modulated referencebeam and phase-sensitive detection of the interference signal.
 14. Amethod in accordance with claim 1, further comprising interposing areference reflector along the local propagation axis.
 15. A method inaccordance with claim 1, wherein the step of interposing a referencereflector includes interposing a microscope coverslip.
 16. A computerprogram product for use on a computer system for determining thethree-dimensional susceptibility of a sample, the computer programproduct comprising a computer usable medium having computer readableprogram code thereon, the computer readable program code including: a.an input for receiving an interference signal obtained throughsuperposition of scattered light from the sample with a reference beam;b. program code for deriving a forward scattering model relatingmeasured data to structure of an object; and c. program code for solvingan inverse scattering problem based upon the forward scattering modeland the interference signal to infer a three-dimensional structure ofthe sample.